{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <div align='center'> 协方差矩阵 </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import matplotlib.pylab as plt\n",
    "\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm\n",
    "\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "plt.style.use('ggplot')\n",
    "plt.rcParams['figure.figsize'] = (8, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 基础"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## np.cov\n",
    "\n",
    "def cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None,aweights=None)\n",
    "\n",
    "- m: 默认每一列代表一个观测点, 行代表属性特征, 由rowvar控制\n",
    "\n",
    "- rowvar: 默认True, 表示行是变量(属性特征)\n",
    "\n",
    "- bias: 默认Flase, N - 1. 这就是为啥np.var()的值和np.var对角线上的值不同的元凶.\n",
    "\n",
    "\n",
    "协方差矩阵计算的是不同维度之间的协方差，而不是不同样本之间的\n",
    "\n",
    "拿到一个样本矩阵，首先要明确的就是行代表什么，列代表什么"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## np.cov(X) 和 np.cov(X, Y)区别"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 5, 6],\n",
       "       [4, 3, 9],\n",
       "       [4, 2, 9],\n",
       "       [4, 7, 2]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.array([[1 ,5 ,6] ,[4 ,3 ,9 ],[ 4 ,2 ,9],[ 4 ,7 ,2]])\n",
    "X\n",
    "# 默认: 行是特征值, 列是观测值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = X[0:2]\n",
    "y = X[2:4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 7.        ,  4.5       ,  4.        , -0.5       ],\n",
       "       [ 4.5       , 10.33333333, 11.5       , -7.16666667],\n",
       "       [ 4.        , 11.5       , 13.        , -8.5       ],\n",
       "       [-0.5       , -7.16666667, -8.5       ,  6.33333333]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.cov(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 7.        ,  4.5       ,  4.        , -0.5       ],\n",
       "       [ 4.5       , 10.33333333, 11.5       , -7.16666667],\n",
       "       [ 4.        , 11.5       , 13.        , -8.5       ],\n",
       "       [-0.5       , -7.16666667, -8.5       ,  6.33333333]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.cov(np.vstack((x, y)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 7.        ,  4.5       ,  4.        , -0.5       ],\n",
       "       [ 4.5       , 10.33333333, 11.5       , -7.16666667],\n",
       "       [ 4.        , 11.5       , 13.        , -8.5       ],\n",
       "       [-0.5       , -7.16666667, -8.5       ,  6.33333333]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# X = np.vstack((x,y))\n",
    "np.cov(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-3.,  1.,  2.])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(np.power(X[0] - np.mean(X[0]), 2)) / (3 - 1)\n",
    "\n",
    "(X[0] - np.mean(X[0])).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实现"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cov_matrix(X, bias=False):\n",
    "    # R是属性特征, C是观测点\n",
    "    R, C = X.shape\n",
    "    if bias:\n",
    "        dof = C\n",
    "    else:\n",
    "        dof = C - 1\n",
    "    cov_matrix_result = np.zeros(R*C).reshape(R, C)\n",
    "    for r in range(0, R):\n",
    "        for c in range(0, C):\n",
    "            cov_matrix_result[r][c] = np.sum((X[r] - np.average(X[r])) * (X[c] - np.average(X[c]))) / dof\n",
    "    return cov_matrix_result "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 7.        ,  4.5       ,  4.        ],\n",
       "       [ 4.5       , 10.33333333, 11.5       ],\n",
       "       [ 4.        , 11.5       , 13.        ],\n",
       "       [-0.5       , -7.16666667, -8.5       ]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cov_matrix(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 残差协方差矩阵"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 模拟数据, x[1,3,5,7,9] 5个值的150观测点, 一个X下10个样本点\n",
    "\n",
    "# np.random.seed(520)\n",
    "# x = np.random.uniform(0, 10, 30)\n",
    "\n",
    "x = np.zeros(5 * 10)\n",
    "x[0:10] = 1\n",
    "x[10:20] = 3\n",
    "x[20:30] = 5\n",
    "x[30:40] = 7\n",
    "x[40:50] = 9\n",
    "\n",
    "e1 = np.random.randn(50)\n",
    "y1 = 2*x + e1\n",
    "\n",
    "e2 = e1*(x + 5)\n",
    "y2 = 2*x + e2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f66f2ae6320>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD7CAYAAABkO19ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHr1JREFUeJzt3X1wHPd93/H37h4Oh+OBJGDwwaAEShZFuaSmpkVRmtKmFDmdOCPHltNKv44djSLVNj2KXdepXU/tdOSH1G3GsmesNp6kjJUqji2nv6i13E4VR5P6QRpTGlOKQDtUJZG0KJCmSIICH3AEcA+72z/2wAcQJJe4A3Zv7/Oa4Rz4E7D71RH44rff35MThiEiIpJdbtIBiIjI/FKiFxHJOCV6EZGMU6IXEck4JXoRkYxTohcRyTglehGRjFOiFxHJOCV6EZGMy7XiIsaYjwNfB9YAy4GvAh7wjLX20zEuoeW5IiKXz4n1Sc1ugWCMuRr4U6AHuAd4HHiftXbEGPN3wAPW2u2XuEx48ODBpuIAGBgY4OjRo01fp5UUU3xpjEsxxZfGuLIc0+DgICxEojfGOMD/Bv4NsA34/cbru4EHgfXAd621D83ytVuBrQDW2o3VanXOcUzL5XLU6/Wmr9NKiim+NMalmOJLY1xZjimfz0PMRN9s6WYr8ENr7SvGmOm2QeDbwGeBW7jAOIC1dhvRLwWAsBW/4bL827uV0hgTpDMuxRRfGuPKckyNHn0szSb69wKLjTHvAzYAf0KU6DcBrwN/BDzQ5D1ERKQJTSV6a+1vTX9sjPkxcC/wFuB7QA140lq7o5l7iIhIc1oy6wbAWvtrjQ/3ATe36roiItIczaMXEcm4lvXoRUTkIqpVvJEROHAAb2oKf2gIopkz8049ehGR+Vat0jU8jFsu4wQBbrlM1/AwtGBaeRxK9CIi88wbGYl6724j5bou5PNR+wJQohcRmWfO1NSZJD/NdaP2BaBELyIyz8JCAYLg3MYgiNoXgBK9iMg884eGonr8dLIPAqhWo/YFoFk3IiLzLZ+ntmED3sgIoesSlEoLOutGiV5EZCHk8/hr1sDAAP4C77+j0o2ISMYp0YuIZJxKNyLSnARXfEo86tGLyNwlvOJT4lGiF5E5S3rFp8SjRC8ic5b0ik+JR4leROYs6RWfEo8SvYjMWdIrPiUezboRkblLeMWnxKNELyLNSXDFp8Sj0o2ISMYp0YuIZJwSvYhIxinRi4hknBK9iEjGKdGLiGScEr2ISMYp0YuIZJwSvYhIxinRi4hknBK9iEjGKdGLiGScEr2ISMYp0YuIZJwSvYhIxinRi4hkXNMHjxhj/hbIAwVgK9ADfBXwgGestZ9u9h4iIjJ3TfforbXvttbeBvwNsBnYBtxtrX0HsMEYs7nZe4iIyNw5YRg2dQFjzHuALwJLgVuBx4F3Aw8C64HvWmsfmuXrthI9AWCt3VitVpuKAyCXy1Gv15u+TisppvjSGJdiii+NcWU5pnx0Lq8T53ObTvTTjDG3A/cDNwA7gc8CtwDubIl+hvDgwYNNxzAwMMDRlJ1ZqZjiS2Nciim+NMaV5ZgGBwchZqJv5WDscWACKAMfBn4O3A5sb+E9RETkMjU1GGuMuRp4BJgCTgKfBK4DvgfUgCettTuajFFERJrQVKK31r5KVJc/26+Am5u5roiItE7T0ytFpMNVq3gjI3DgAN7UFP7QEEQDhZISWjAlInNXrdI1PIxbLuMEAW65TNfwMLRgFp20jhK9iMyZNzIS9d7dRipxXcjno3ZJDSV6EZkzZ2rqTJKf5rpRu6SGEr2IzFlYKEAQnNsYBFG7pIYSvYjMmT80FNXjp5N9EEC1GrVLamjWjYjMXT5PbcMGvJERQtclKJU06yaFlOhFpDn5PP6aNTAwgJ+y7QYkotKNiEjGKdGLiGScEr2ISMYp0YuIZJwSvYhIxinRi4hknBK9iEjGKdGLiGScEr2ISMYp0YuIZJwSvYhIxinRi4hknDY16yQ62zM+vVeSIerRdwqd7Rmf3ivJGCX6DqGzPePTe3WZqlW8PXtgeDh61S/E1FGi7xA62zM+vVeXQU8/bUGJvkPobM/4wkIBKhXcAwfg5Zej10pF79Us9PTTHpToO4TO9ozPX7mS3K5duOPjEIa44+Pkdu3CX7ky6dBSR08/7UGJvlM0zvYMSqXTZ3vWNmzQTJJZeIcOUV+/nqBUAschKJWor1+Pd+hQ0qGljp4U24OmV3YSne0ZizM1Bd3dBFdeCUuXEhw/fqZdzuEPDeEOD5/pMOhJMZXUoxeZQb3Uy6AnxbagHr3IDOqlXiY9KaaeEr3ITI1eqjcycrqXqpWxF6FVxKmnRC8yG/VS42nMoyefx8nncctl3OFhlW9SRjV6EZkzzaNvD0r0IjJnmkffHlS6EZE5CwsFnHL53GQfBITFYnJBpVWCYxnq0YvInGnFdUwJ7wnUVI/eGFMEtgGrgCJwD7AE+CrgAc9Yaz/dbJAiklKaoRTLxcYy/DVr5v3+TfXorbUTwOettbcBDwO/R5T477bWvgPYYIzZ3HyYIgtMW+/GNz1DacOG6FVJ/jxJj2U0XaO31u5tfDgIjAE1oGyMeRgoAZuA7TO/zhizFdjauAYDAwPNhkIul2vJdVpJMcWXmriqVdixA6e7G89x6M/lCF99FTZtSkUSS837NEMa40pNTIODnBg5yd/90OP4CZelS/r4p+/yWXrFYliA+JwwDJu+iDHm/cB9wIeBYWAn8FngFsC11j50iUuEBw8ebDqOgYEBjqZszrNiii8tcXl79uAeO4Y7Osriri5O1moEy5YR9PUtyGP2paTlfZopjXGlJaZjh2vYf/sKdbebQrHA1MQUuaCCeXAtfSu65nTNwcFBACfO5zY9GGuM+RBwF3CXtXYUKBMl/J8DtzNLb14kzZzxcXK7d+OeOhVtU3zqFLndu3HGx5MOTdrU408sYW//jVS6ewlwqXT3srf/Rh5/YsmC3L/ZwdgbiWryPwWeNMZUgfuB7xGVcJ601u5oOkqRBeQeOwaed+7AmedF7SJz8MYbLuQ9DuXXcrxQYKpRm3/jDX9B7t9UorfWPkc0u2amm5u5rkiSgr4+vMOHz50yGAQEfX3JBiZt601vCnj9dQ/vrGzp+1H7QtCCKZEZwt5e6mvW4B45Eh08UiwSLF9O2NubdGjSpu64Y5KfPxcy9coB6vUqlVye/NoruOOONpl1I5I1/tAQ7tgYwapV0N9PMDamRUDSFKdWZfmeFzj6yjieX8Xx8gy4+3Fq64G5DcZeDq2MFZlJh2lIi/31l19n6h8O0O1P4rnQ7U8y9Q8H+Osvv74g91ePvpNo3/D4tE2xtNDeJw/gOXnCxgB/6LrUgjx7nzwAzP+Tonr0nSLhvTZEOlm95hDiEAYOgQ9hEP29Xos1Db5pSvQdQvuGXyZtgSAtdHzganJhhSAMCIAgDMiFFY4PXL0g91ei7xBJ77XRVvT0Iy22Ystq9rCGCYoEuExQZA9rWLFl9YLcXzX6DhEWCjiNZf0cPozbWNYfam74eZLeaVCyZ+/+RQzn/wlvru6jyCQT9PB6/iqC/Q5Qmff7q0ffIfyVK8nt2oU7Ph4t6x8fJ7drF/7KlUmHljp6+pFWO3jQxfcdXBxcwMXB9x0OHlyYFKxE3yG8Q4eor19PUCpFi4BKJerr1+MdOpR0aKkTFgpnVsVOC4KoXc6n8YxLcmo13u7voJeTeE5ALyd5u78Dp1ZbkPurdNMhnKkp6O4muPJKWLqU4PjxM+1yDn9oCHd4+MzUU52adGHVKlNP7+RH23uZqoQUuuvctnknhS1v09Tds6wOXuWUmycMHBwgxMF386wOXgX65/3+SvQd4nSN/siRqEZfrUbL+lWjP59OTYqt/Iv9fOvbfQS49PQ4TE7mefWXfdyzdD+ljdckHV5qDC0rs/sNF8cPCUMHxwnxPIehZWUWItGrdNMhTtfoy+WoRl8uq0Z/MTo1KZanfhAS4J4zbh3g8tQPmj/nIkuu3+TRt8SnWIRFi6BYhL4lPtdvmm1PyNZTj75DTNfo3dHRqEbf20vwlrfgHTqkmSQyZ6OnilA5wdT+Ufx6jVqui/yVyxg9tTD7rLeL37x/OUdfeoX9R3oIwjyuU+XK5ZP85v1rF+T+SvQd4vRMkukTxcJQM0mkaZWBN8P3tzNRLRG4XbjBBPmxXVRu/OdJh5YqfSu6uPuhtTz1F0fwx8Hr7eaW371yzqdLXS4l+g4Reh5dL70UlSC6unAnJnBfeonqDTckHZq0sfzoQX488XaW+4foDqtUnF5erF/Lb4weBFYlHV6q9K3o4o7PrErkeEMlepHZaAO4WF4ZruN73ewPr4LGfJKcF/LKcD3hyORsGoztEI7vU3/rWwkWLYpq9IsWUX/rW3H8hTnKrK1oC4TYDh7vJQwCgiA6MSkIIAwCDh7XIS1pokTfIcJCATyP4Ior4LrrolfP0yKgWWgDuPiO962GSo3ADwlDCPwQKrWoXVJDib5D+ENDUY/07HNQtQhoVtoCIb7xSp7n3U2Msxg/dBlncfT3ispcaaIafafQIqDYwkIBp1w+N9kHAWGxmFxQKdXTA+S7eHlqLdM1+kI+pKdHNfo0UaLvJDo1KRZ/aAh3xw7cY8eiVcSVCkFfH/66dUmHljqFgs/UVBdRkgdwmJqK2iU9VLoRuZCz1xzIrA4dOjvJT3Ma7ZIW6tF3Ek0ZjMUbGYFFiwh6e89sABcE2o9+FmNjLrlciO+fSfaeFzI2pj5kmuhfo1NoymBsGoyNr6sretrJ5UK6uqLXs9slHZToO4Q3MgKOg3vwILz8cvTqOJoyOAvtRx/fb//2JHB+lWu6XdJBib5DOOPj5Hbvxj11Ktq98tQpcrt344yPJx1a6mgqanyf+ESZW2+dpFQK6OoKKZUCbr11kk98opx0aHKWbNToVXu+JPfYMfC8cxcBeV7ULufSVNTY+vtDvv71k3z/+z1MTJQoFsvcccck/f0q3aRJ+yf6Ru2ZfB4nn8ctl3GHh6lt2KAfzLMEfX14hw+f20sNAgIdPDI7TUWNrb8/5L77JhgYKHL06ETS4cgs2r50o+Xq8YS9vdTXrCEoFqO9bopF6mvWEPZqT5LZ7N3r8tGPLuXXf93jox9dyt69bf+jIh2s7b97NUMiHn9oCMKQYNWqaK+bVasgDFV3nsXevS733vsmXnyxi9FRhxdf7OLee9+kZC9tq+2/czVDIqZG3TkolU7XnVXemt1XvrKYqSk4fNjltdei16mpqF2kHbV9jd4fGsJt1OgBzZC4GNWdY9m3z+PIEQ/fd3AcCEOHSsVh377g0l8skkJtn+g1Q0JabXTUpV6PVnpOzwuv1x1GR9v+AVg6VPsnelBPVVoqCGbu3XLxdpG0UxdFZIaJidkT+oXaRdKuqR69MaYAfAT4FPAFa+0jxpibgK8CHvCMtfbTzYcpsnB6ewMmJhzC8Exid5yQ3l7V6KU9NdujXwFMAt85q20bcLe19h3ABmPM5ibvcWnVKt6ePTA8HL1qoy5pwnXX1Ru1+enVndExedddp8M0pD05YQv22jbGfAHYB3wfeBJ4N/AgsB74rrX2oVm+ZiuwFcBau7E61+RcrcKOHTjd3XhdXfi1GmGlAps2pWJANpfLUa+nK0GkMSZIT1yf+YzLN74RDchGs26iXRk/9rGAr3wl+V59Wt6nmdIYV5Zjykf5LVY9cT4GYweBbwOfBW7hAk8N1tptRL1/gPDoHAdRvT17cE+exB0dZXFXFydrNYJlywj+/u9TsXf4wMAAc/1/my9pjAnSE9dzz/WTz0cHZ4Shg+OE5PMhzz1X5+jRsYSjS8/7NFMa48pyTIODg7E/t6WDsdbaY0AZ+DDwc+B2YHsr7zGTdmWUVjtyxMXzoFgMKZWiV8+L2kXaUbODsVcAjxP14ivGmHcB9wPfA2rAk9baHU1HeRHalfEyaJfPWFasCDh0yJu1XWah76vUayrRW2sPADfO8p9ubua6l0O7MsZUrTL19E5+tL2XqUpIobvObZt3UtjyNv1QzrBlS4X9+z3K5TM9+FIpYMuWSoJRpZR2j20Lbb9ganpXRvfIkdO7MgbLl2tXxhnKv9jPt77dR4BLT4/D5GSeV3/Zxz1L91PaeE3S4aXKBz84wUsv5XjtNY8g6MJ1a6xe7fPBD2oL3pkutntsGsbIJNL2id4fGsIdG4t2Y+zvJxgb0143s3jqByEB7jk/j0Hg8tQPAm7fmGxsadPfH/KlL00fppGjWJzSYRoXoN1jL0OCJa62T/Ta6yae0VNFPGeC8Kzxd88JGD1VTDCq9OovVfjwlt30FQocm5rCLw0B+p6aKSwUcI4di56oDx/GrVajJ2qVTs+VcIkrG9MIpve62bAhelWSP9/VQ0yeqDD8gsfTTzsMv+AxeaICV+vJ5zyNH0q3XMYJAtxyOfoh1UK88/grV5LbtQu3XI5mvZXL5Hbtwl+5MunQUiXpA5Kykejlkq6/IeSbO9/J3iNLeeOYx94jS/nmzndy/Q0qR8yU9A9lO/EOHaK+di2UyzAyAuUy9bVr8Q4dSjq0VEm6xNX+pRuJ5WtfW0y52s1L4VtPt3nVkK99bTGPPqqpqGdL+oeynTjj4+T27YPFi2HJEjhxgty+fdR6epIOLVXCQgGnXD73+yoICIsLUzrNRo9ee91c0gsv5PH9c1dL+77DCy+ozDWTTi2LT+tY4vGHhqK8dPY08AWcNNL+iV711FgmLjAz8ELtnSzpH8p2EvT1nV67EjVoHcusEj7Ks+1LN5rHG0+hEFIuz94uM2gmV2xax3IZEjwgqe0Tveqp8Sxf7jdWep5dvglZvtxPKqR006llsWgdS3to+9KN6qnxbNlSpacnxPOCRhk1oKcnZMsWlbikCfk8tXXrcMbGCHfvxhkbo7ZunZ5+UqbtE73qqfFce63PbbdVWLUqYPnykFWrAm67rcK116pHL02oVul68UXC/n6ca68l7O+n68UXNUaWMm1fulE9NZ477phk164uNm+usmhRgVOnqjhO1C4yVxojaw/tn+hB9dQY+vtDPve56f1buikWa9q/RZqmMbL2kI1EL7H094fcd98EAwNFjh7VvEppXtILgSSebCR6HXwQy9iY0+jRuxSLRfXopWn+0BBuY7MuQGNkKdX+iV4HH8QyNubwwAOLG3usO7hugeef7+JLXzqpZC9zpzGyttD2s260AVU8jz5aZOfOPCdOeExOOpw44bFzZ55HH9Uj9qy0rUZ82j029do+0WswKJ5nn83jOOA01ktNf/zss/qhPI+21ZCMafvSTVgocPJXp/jJ0z1MTTkUCgVu3TJJ7yotmDpbeIHqzIXaO5mmDErWtH2P/mjvav7y4Ry/3OMyNubwyz0uf/lwjqO9q5MOLVU2b64ShmcS+/THmzerlzqTnhIla9q+R//4E0sY7trEiZ0HyFUq1Lv7WfK2K8g/AffdpymE0z7wgejA6337cgSBh+sGXHVVnQ98QO/RTJoyKFnT9ol+926P//t0kXr9HxGGLs6pgNzT0LdCCexs/f0hX/zi9IIpj2JxUtMrL0BTBmVe6HDwuXvmmTyVSjTC6LrRYRq+H7XLubRgKiZNGZRWS3gaeNsn+okJ93S9eXpfM8eJ2kXmTNtqSAslPcDf9om+qyskT5WhYB9FJpigyIh3FV1dSvQikg5JD/C3fTa85soJNgY7WMxJPCdgMSfZGOzgmitVmhCRdEj63Iy2T/Q39O8lzHeB60SLgVyHMN/FDf17kw5NRARI/tyMti/dTL5RYVHJo34iJAwdHCdkUclh8o1K0qFJO9NGedJKCQ/wt32P3lnUzYljNJJ89HriWNQuMifaAkEypu0T/U9eW0uPVyH0A+p1CP2AHq/CT15bm3Ro0qa0UZ60XMKdh7ZP9DUnz/bazZxgCQEuJ1jC9trN1Bw9ZsvcJD1DQrIn6c5D29fox8ZcfM9jL9cRhtO7M4aMjenQa5kbbYEgrZZ056Hte/TLlgW4LjhOSC4Xvbpu1C4yF0nPkJDs0fTKJq1e7bN8uU93d4jnQXd3yPLlPqtXq0cvc9SYIRGUSqdnSOjEMmlG0p2HlpdujDEe8J+B6wEPuN9a+4tW32faZz5zknvvfRMrVgR0dTnUasHpdpE50xYI0koZnF55J+Baa28F/h3wtXm4x2nXXBPwyCNvsG5djWXLQtatq/HII29wzTUq3YhIiiR45KITtviIIWPMQ8CTNHrzwDpr7XmngBhjtgJbAay1G6stmGaUy+Wo1+tNX6eVFFN8aYxLMcWXxriyHFM++kXhxLpn03eb3b8H/gZ4H/DSbJ9grd0GbGv8NTzazONxuUz+Zz9jcRhy0nGo3nQTlEpzv14LDQwM0NT/2zxIY0yQzrgUU3xpjCvLMQ0ODsb+3Pko3TwPvGCt/RJwI/D/5uEeZ5TL9Dz2GN7oKI7v442O0vPYY1Auz+ttRUTaxXwk+u8Ci40xPwH+A/D783CP0/I/+xn09ECu8XCSy0FPT9QuIiKtL91Ya2vA3a2+7oW4J0+eSfLTcrmoXURE2n8efbB4Mcwc2KjXo3YREWn/RF+96SaYnDyT7Ot1mJyM2kVEpP0TPaUSk3feib9sGaHn4S9bxuSdd6Zm1o2ISNLaflMzAEolqu96FwwMUE3ZVCoRkaS1f49eREQuSoleRCTjlOhFRDIuGzV6HeQsImmXYJ5q/x69DnIWkbTTmbHN8UZGwHFwf/UrePnl6NVxdJCziKRG0mfGtn2id8bHye3ZgzsxAWGIOzFBbs8enPHxpEMTEQF0ZmzT3GPHojfw7N+Urhu1i4ikgM6MbVLQ1we+f+5ZjL4ftYuIpEDmzoxdaGFvL/Vrr8UdHQXHIVi0iGDZMsLe3qRDExGJJHxmbNsnen9oCHdsjGBwEPr7CcbGFvQ3pYhILAkeON/2iT7p35QiImnX/okeEv1NKSKSdtlI9CIiaaeVsSIiGaaVsbJgqlW8PXtgeDh61TYRIgtCK2NlYWhPIJHEaGWsLIikexQinUwrY2VBJN2jEOlkSa+MVaLvEEn3KEQ6WmO9T1AqnV7vU9uwQStjpbX8oSHc4eEz31gL3KMQ6XhaGSvzTiuIRTqWEn0n0QpikY6kRN9JdLauSEfSYGyn0Dx6kY6lRN8hNI9epHMp0XcIzaMX6VxK9B1C8+hFOpcSfYdIemWeiCRHs246hebRi3QsJfpOonn0Ih1JpRsRkYybc4/eGHML8Dmg21p7W6NtEfBnwCDgA79rrT3QikBFRGRumunR/2PgDwHnrLaPAzuttb8G/DfggSauLyIiLXDJHr0x5r8CG2c0f8ha+8fGmKtmtG8G/sAY8yHg/cAVF7nuVmArgLWWgYGBy4l7VrlcriXXaSXFFF8a41JM8aUxLsXUuOelPsFa+9HLvOYfA/8d+B3gJxe57jZgW+Ov4dEWDA4ODAzQiuu0kmKKL41xKab40hhXlmMaHByM/bmtHox9HnjCWvsnwHuAn7b4+iIicpmcMAzn9IXGmG8CNwFXAy8D/wI4Cnwb6AXKwL+01h6Jcbm5BSEi0tmcS38KEIZhZv7cddddzyUdg2LKVlyKqb3jUkzRH82jFxHJOCV6EZGMy1qi33bpT1lwiim+NMalmOJLY1yKiSYGY0VEpD1krUcvIiIzKNGLiGRcJrYpNsYUgI8AnwK+YK19JNmIwBhTJKrFrQKKwD3W2peTjQqMMX8L5IECsNVa+4uEQwLAGPNx4OvAGmvtvoTDAcAYMwU82/jrQ9ba7yUZD4Ax5ibgi0A38B1r7cMJx3Mn0R5XEOWTFdbaaxMMCTidE75FtMFiAfiP1tr/mXBMXcDDwOpG09aFygmZSPTACmAS+E7SgUyz1k4YYz5vrd3b2Nfn94B/nYK43g1gjHmAaG+ixBO9MeZq4L3A9qRjmeFQY4O+VDDG9AIPAe+z1o4mHQ+AtfYx4DEAY8xHgDcnG9FpNwMVa+07jTEbgc8DiSZ6om1h9ltr72ns/vuHgFmIG2ci0VtrXwO+aYz5QtKxnM1au7fx4SDwyyRjmWaMeQ9Rj3ApcEvC4WCMcYD/Avwr0jdDom6MeYpoxfenrbVJ/xtuIXoa+4tG0n/QWvu/Eo4JON2D/jhRjGnwLPDJxgr+1cAfJBwPQA9nVrLuBq5fqBurRj/PjDHvB94OfCPpWACstf/HWnsj8AngT5OOh2gH0x9aa19JOpCZrLVrrLW3EG3r8eWk4wEWA89Za28H/hlR7z4tPgb8lbX2ZNKBNCwGFgE/A2rApmTDAaJS0pAx5kdE/3anFurGSvTzqLFd813AXdbaetLxzHCcqNyVtPcC7zfG/BjYAPyVMea6ZEM6jwMcSzoIol7gdH23ygImiosxxiwGPkz0ZJYW7wN+3tgldyvRGF6irLWnrLV3Nw5q2s4ClpozUboxxlwBPE5UIqkYY95lrb0n4ZhuJCpF/BR40hhTtdb+RsIxXQ08AkwBJ4FPJhkPgLX2t6Y/biT7e9MwGGuMWQ78D6ACHCHqsSbKWvu8MeYlY8x2ol5q4mM+DZ8C/txaW046kLM8AfxO43sqR1QPT5QxZiVRcneAHwL/aaHurQVTIiIZp9KNiEjGKdGLiGScEr2ISMYp0YuIZJwSvYhIxinRi4hknBK9iEjGKdGLiGTc/wdLZhPjexjprAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# sns.regplot(x, y1, color='b')\n",
    "# sns.regplot(x, y2, color='r')\n",
    "\n",
    "plt.scatter(x, y1, color='b', alpha=0.6)\n",
    "plt.scatter(x, y2, color='r', alpha=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "homo_fitted_result = sm.OLS(endog=y1, exog=x).fit()\n",
    "hetero_fitted_result = sm.OLS(endog=y2, exog=x).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Residuals vs. Fits Plot')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxsAAAF0CAYAAABL4yRvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3X98pFd92PvP49Fgs1hE1k6QvJARu9AQWqDQkEA2/Myr1/eGS2g23R7U2MUkEBNAuZvLkmmccMtCkm47r5h2iVKKA4VtcBGHpCYNNYU2KaRN1jRB1CEpSWnWzNSsJSMrG8ssNtL4uX/MyJbk0e5IOzPPo5nP+/XSa3fOPDPPd2dHc+b7nHO+J0nTFEmSJEnqtiuyDkCSJEnSYDLZkCRJktQTJhuSJEmSesJkQ5IkSVJPmGxIkiRJ6gmTDUmSJEk9YbKhnkiS5ElJktydJMlsB8d+NkmSE1067wuTJEk3/Dy9G8+7lyRJciRJkvNJkrx0m/tHkyQpZhDX7Ib/l8/u8jmuSJJkrMuhSdoj7FuyY9+i3TLZGEKtD+A0SZJHkiS5N0mSDyZJMtrl06wBdWCxy897UWma/lGapglwsJ/n7Zctnd36z3dtOeyvgP8NPNjm8U8D/hJ4ag9i++yGmO5LkuQ3kiR55vr9aZrOtP5v3nUZp/kU8NOXHaykrrNv2bvsW+xbeslkY3i9CygArwC+D3hnN588TdOH0zR9eZqmv9DN5xUAr0zTNNnw82cb70zT9HfTNH1umqZfbPPYEZr/773yrtaH/ncDDeA/JEnyhC4+/5VdfC5J3WffsnfZt6gnTDaGWNr05zQz+kevYCRJUkyS5J8mSXJ/kiRLSZLcliRJacP91yRJ8q9bw6mrSZJ8LUmS79xw/yc3XIX48NbzJkny9CRJPpEkyTeSJFkBvnfDfa9IkiTdcvxXkyR5fevvT0yS5D8kSfL1JEkeTpLkfyRJ8iM7+XcnSfKPkiQ5lyTJWuvfeEOHjwutc16zoe1g6yre97Rul5Mk+a0kSR5MkuRbSZLUkiS5eifx7VaSJKUtV6VeseG+p7de17tbTXdvc9ybWlMUVlv/vz+7m1jSNP3fwAzwDOBvdhj/05Mk+e0kSS4kSXJPkiS/miTJk1v3fbYV/8uBd17s/SUpW/Yt9i32LdrIZGOIJUlSSJLk+cCPAP9xw10ngZcAfwt4HnAN8Gsb7n8nUAaeDTwJeBlwbv3ONE1f3boCcbrNOa8Cfgf4MvA0mh8Yf76DsFPgw624xlp//zedfugmSfLDwE8Bfxt4IvCcVjyd+ATwADC9oe11wF1pmv5h6/Z7gQs0X59R4AeBhzp8/suSpunS+hWpNvd9dcsUgIMbrl59FiBJkhcAs8CP0vx//WvAv7mMkEZaf379Uge23he/S/O98B3A/0nzPfjeVvyvaMX/OVpXuFo/r7+M+CT1gH2LfYt9izYy2Rhe76Q59/UPgVNpmp6C5iIp4M3Az6ZpWkvT9Fzr2B9KHhuyfJDmvMxnA2tpmv5FmqaPm8O5jR+h+UH882ma/mWapvfRnAfakTRNH0rTdC5N03vTNP0m8EGaw5+HOnyKB4Grge8Biq3nubfDc38L+NfAjQBJkiQ0O4Rbtzz/M4BntIb7/0eapmsdxtap/7zh6kunr3snLrT+fDGwL03Tr6dpWt/pkyRNT6fZudyRpulXO3jYq2m+L/5hmqb3p2n6p8DNwHTrPSlpb7BvsW/Zyr5lyPlCD693ASXgszTn1a4rAfuA31v/0AH+G825mJOtY34RmAM+AnwtSZJ/nCTJvg7P+3Tgq2maPrKboJMkuTJJkl9pDYWuAUutuzp6L6dp+p+AnwSOA19PkuTXkySZ2kEIHwBelCTJs4CX0nxNbttw/wzweeDTSZKcTZLkZ5Ik6fY81o3zars2jN6a9vDDQAAWW9MRnrvDp3knzfm0d9JcwPmjHT7uEM33RWND21manf3+HcYgKTv2LfYtm9i3yGRjiKVpej/wFuBokiQvaTXfD6wC37dloViyfiWidQXo/6N5BepHgb9P51UgltjwC54kyZVsrl6x1movbPjzqg33vw34v2guPizSHIbfkTRNT6dp+jyaC80mgN/YwWO/DPw+zatOrwM+lqbpAxvuP5+m6U8BT6FZ2eLnaHZAmyRJcqAHHUUn1j9wR9rdmabpv0/T9PuAZ9K8kvafWlfZOvWuNE2vSNN0Mk3Tt6Rp2umVxXuBqS1Xmr4TWEnTdONQeWO72CXlg32LfctW9i3DzWRjyKVp+hXg48A7WrcbQAR+OUmSv5EkyUiSJN+WJMlfW39MkiQvSJJkguYVqS/TnDfZaVWI3wGemSTJ32k9R2Tz1YWzwCPAS1pD67fS/HBdV6JZXu8+mgsP/+VO/r1JkjyztVjsCcDXgP+1g9jX/RrNTvDvsnmYmyRJXpQkyXjr5pdoDuM/Ycsx061zf3iH5+2GczQ/6P9BkiRXJUkyvj4nOUmSpyVJ8p2tOa5LNP9vi/Tnc+JTNDv+k63329+geZVzay39/wm8qtWhXpkkyeTWJ5KUPfsW+xb7Fj0qTVN/huyH5vD2iQ23n0fzQ/h7Wre/DXgfzQ+FVZof+L+y4fj3A8s0rxTdD3wIGNtw/3mai+02/vynDff/ZOvxCzSv0Hx4Szw/S/MD/yzNjur3gde37isD/xX4JvAVmleAzgPPb93/023OnQIvad3/Bpofit8CvkFz8eLf2OHrt691zj9uc9+/p9kJrNG8ovIe4Motx1xHczHgP9nF/10KvGKb+354m3/7DVuO+7s0P1i/1fr/W/9//0Ga9esfprnw8A+Al+72fdXm/v/eJrb/teH+F7X+bx+kOUz+bppznzc+x1OBz7SO+SZQzfr3yR9//Gn+bP0MwL7FviW1b/EnJWm9yJIkSZLUVU6jkiRJktQTJhuSJEmSesJkQ5IkSVJPmGxIkiRJ6olBrCnsindJ2mwn9ezVnn2LJG3WUd8yiMkG586d29XjSqUSS0tLlz6wj/IWU97iAWPqVN5iyls8MJgxHThwoIvRDLdB6lsgn3HlMSbIZ1zG1Lk8xrXXY9pJ3+I0KkmSJEk9YbIhSZIkqSdMNiRJkiT1hMmGJEmSpJ4w2ZAkSZLUEyYbkiRJknrCZEOSJElST5hsSJIkSeoJkw1JkiRJPWGyIUmSJKknRrIOQHtHvV7g+PECtdp+JicbVCorlMuNrMOSJEkCoFCvM1qtUlhYoDE5yUqlAqVS1mENNZMNdaReLzA9PU6tVgAKAMzPF5mbWzbhkCRJmSvU64xPT1Os1R5tK87Pk3760zA6mmFkw81pVOpItTpKrVbc1FarFalW/eWVJEnZG61WNyUaAMVajcKJE9kEJMBkQx1aWCi0bV9cbN8uSZLUT4WFhbbtyb339jkSbWSyoY5MTrafKjUx4RQqSZKUvcbkZNv29Npr+xyJNjLZUEcqlRWmplY3tU1NrVKprGQUkSRJ0mNWKhVWp6Y2ta1OTdFwGlWmXCCujpTLDebmljl1qkS9vsbEhNWoJElSfjTKZZbn5prVqBYXaUxMsFKpcM3Bg7C0lHV4Q8tkQx0rlxucPt1gaen+rEORJEl6nEa5zPnZ2azD0AZOo5IkSZLUEyYbkiRJknrCaVSSJEnSkGi3y3qjXO7Z+Uw2JEkDLYSwD7gVeCqwD3gd8G3ALwMF4EyM8e3ZRShJ/bHdLuvLc3NQKvXknE6jkiQNtBjjBeCdMcZXAh8E3kIz+bghxvj9wPNDCIezjFGS+mG7XdZHq9WendORDUnSwIsx/kXrrweAZWAVeDCE8EHgauB7gD/Y+JgQwk3ATa3HU9rlVb+RkZFdP7aX8hhXHmOCfMZlTJ3LY1xZxTSyvNy2/arlZehRTJkkG+2GtGOMf97muALwXuA5NIe63xxj/FI/Y5UkDYYQwg8DLwDeSDOJ+AhwM/Ay2oz0xxhvpdlXAaRLu6zTXyqV2O1jeymPceUxJshnXMbUuTzGlVVMY+Pj7GvT/tD4OCNrax3HdODAgY7Pmck0qm2GtNs5ClwRY3w58LPALX0KUZI0QEIIbwD+HvD3YoxfBx6kmXT8MfAqtoxqSNIg2m6X9ZVKpWfnzGwa1ZYh7bPbHHYY+GQI4TXAm4FntTtokIe68xZT3uIBY+pU3mLKWzxgTIMqhPBCmiMUvw98JoTwLZp9yu00p1N9Jsb4hxmGKEl9sd0u6wNbjWrDkPbfvchh7wA+BbwG+LN2BwzyUHfeYspbPGBMncpbTHmLBwYzpp0MdQ+qGOMf0ZyKu9WL+h2LJGWt37usZ1aNasuQ9lqr7WgI4eMbDvsC8MUY47uBFwJf7n+kkiRJknYjk2Rjw5D2d9Ac0v5M665J4NkbDv0o8OQQwueAXwT+374GKkmSJGnXMplGtd2QdoxxFpjdcHsVuKGPoUmSJEnqEjf1kyRJktQTJhuSJEmSesJkQ5IkSVJPZFr6VpIkDbdCvd6s+b+wQGFqisKxYz2t+S+pv0w2JElSJgr1OuPT0xRrtWbDmTOMnznD8tycCYc0IJxGJUmSMjFarT6WaLQUazVGq9WMIpLUbSYbkiQpE4WFhfbti4t9jkRSr5hsSJKkTDQmJ9u3T0z0ORJJvWKyIUmSMrFSqbA6NbWpbXVqipVKJaOIJHWbC8QlSVImGuUyy3NzzWpUi4uMlMssW41KGigmG5IkKTONcpnzs7MAlEolGktLGUckqZtMNiRl4u674eabx1hYKDA52aBSWaFcbmQdliRJ6iKTDUl9V68XuP76ImfPPuHRtvn5InNzyyYckiQNEBeIS+q7anWUs2eTTW21WpFqdTSjiCRJUi+YbEjqu4WFQtv2xcX27ZIkaW8y2ZDUd5OT7adKTUw4hUqSpEFisiGp7yqVFQ4dSje1TU2tUqmsZBSRJEnqBReI01ysWq2Osrw8wvj4mFVxpB4rlxvccccqN9+8xuJigYkJq1FJkjSIhj7ZqNcLTE+PU6sVWy37rIoj9cHBgzA7ez7rMCRJUg8N/TSqanV0Q6LRZFUcSZIk6fINfbJhVRxJkiSpN4Y+2bAqjiRJktQbQ59sVCorTE2tbmqzKo4kSZJ0+YZ+gXi53GBubrlVjeoqxscfsiqOJEmS1AVDn2xAM+GYnT1PqVRiacnqOJIkSXtRoV5ntFqlsLBAY3KSlUoFSqWswxpqJhuSJEna8wr1OuPT0xRrtUfbivPzpJ/+NIxaZTQrQ79mQ5IkSXvfaLW6KdEAKNZqFE6cyCYgASYbkiRJGgCFhYW27cm99/Y5Em3kNCpJkqSc27oWgZMnnRq0RWNysm17eu21fY5EG2WSbIQQrgJ+AjgOnIgxfvgixz4E3Nm6eSrGeHvvI5QkScqHdmsR0rvuonDbbTTK5Qwjy5eVSoXi/Pym12l1aorUaVSZympkYwL4JnBbB8cuxBhfcbEDQgg3ATcBxBgp7bLqwMjIyK4f2yt5iylv8YAxdSpvMeUtHjAmSfnUbi1CcvYso9Uq52dnM4oqfxrlMstzc80RoMVFGhMTrFQqXHPwICwtZR3e0Mok2Ygx1oAPhBBOdHD4Wgjh94Al4O0xxrNtnu9W4NbWzXRpl2+oZunbfL0Z8xZT3uIBY+pU3mLKWzwwmDEdOHCgi9FIysJ2axEKi4t9jiT/GuWyCVjO5H6BeIzxmTHGlwEfAX4p63gkSZL6abu1CI2JiT5HIu1crpKNEMLREMLHt7k7Af6yn/Fos3q9wI03Fjh6dD8zM2PU64WsQ5IkaeCtVCqsTk1taksPHWpuWCflXFYLxJ8GfAI4ADwcQviBGOPrgEng2RuOewrwm8DDwH3AWzMIVzQTjenpcWq1AtBMMubni8zNLVMuN7INTpKkAdZuLcLIyZM0rEalPSCrNRv3AC9s0z4LzG64fR/w0j6Gpm1Uq6PUasVNbbVakWp1lNnZ8xlFJUnSkEnTrCOQdsR9NtSRhYX2U6YWF51KJUlSL1n6VntZrtZsKL8mJ9tPlZqYcAqVJEm9dLHSt1LemWyoI5XKClNTq5vapqZWqVRWMopIkqThYOlb7WVOo1JHyuUGc3PLnDpVol5fY2KiQaWy4uJwqYfq9QLHjxeo1fYzOenvnDSsLH2rvcxkQx0rlxucPt1gaen+rEPJNb8gqhusACdp3UqlQnF+fvOaDUvfao8w2ZC6yC+I6hYrwElaZ+lb7WUmG1IX+QVR3WIFOA2LQr3e/BK9sEBhaorCsWNWWGqjUS5zfvbR3QEolUqwtJRhRFJnTDbUMacHXZpfENUtVoDTMHhcSdczZxg/c4bluTkTDmlAWI1KHVmfHjQ3V+DMmSu5/fZ9TE+PU6/7JXojvyCqW6wAp2HQrqRrsVazpKs0QEw21JGLTQ/SY/yCqG5ZrwA3Pd3g8OGHOXLkgmt/NHAs6SoNPqdRqSNOD+qMJYLVTVaA06CzpKs0+Ew2aE4RqlZHWV4eYXx8zC+HbTg9qHN+QZSkzrQr6bo6NWVJV2mADH2y8Vip0vUpQvssVdpGpbLC/Hxx01Qqpwe150J6SepMo1zm/Hvew9ixYxQeeACuuYbz73mPi8OlATL0yYalSjvj9KDOuM+GJHWuUK8z9ra3UbznnmbDAw8w9ra3WY1KGiBDv0DctQg7l6ZZR5BfLqSXpM5ZjUoafEM/suFahM54xb4zJq+S1DmrUUmDb+hHNiqVFQ4c2Fyq9MAB1yJs5RX7zpi8SlpXqNcZm5lh5LrrGJuZoVCvZx1S7liNShp8Qz+yAZAkyUVvyyv2nXIhvZQ/IYSrgJ8AjgMnYowfDiF8L/DLNIdqz8QY397Nc27dGXsfUJyfdy3CFlaj6lyhXme0WqWwsNBM0k6ehFEv+Cn/hn5ko1od5Wtf25xzfe1rI16x38Ir9p1xI7bO3X03zMyMcfTofmZmxtyNXr00AXwTuG1D263ADTHG7weeH0I43M0TuhahM+vVqFaf9jQeefKTeWRqympUbawnr/tuv50rz5xh3+23U3zVqxwt054w9CMbXrHvjFfsO+c+G5dWrxe4/voiZ88+4dE21wCpV2KMNeADIYQTACGEa4BV4MEQwgeBq4HvAf5g4+NCCDcBN7Weg1Kp1PE5R5aX27Zftby8o+fppZGRkexjuftuij/zMyQbqlGVfuZnWL3jDjh4MNvYNsj6tSocP05hS/KanD1L6dQpGqdPZxTV42X9Om0nj3ENU0xDn2x4xb4zlr5VN1Wro5w9u3m6oiWn1WcHgI8ANwMvo81If4zxVpojIADp0tJSx08+Nj7OvjbtD42Pc34Hz9NLpVKJnfybemHs5pt5wtmzm9qSs2dZu/lmzs/OZhTV42X9Wu2v1Wh3CXStXuf+nLyfIPvXCR4/3WylUuGav/W3Mo9rqzy8VlvtJKYDBw50/LxDP42qUllhamrzAnGv2F+cpW91uRxRVJZijH8JPAi8Efhj4FVsGdW4XCuVCqtTU5vaXIvweFaj6owL6TvTbrrZ+PR0c96uMjP0IxvrV+yr1VGWl69ifPwhr9i3YelbdZMjiuqnEMLTgE/QHM14OITwA8CbgdtpTqf6TIzxD7t5zka5zPLcHKPVKlctL/PQ+DgrlYprEbbwS3RnvnHDDTzxt3+bZG3t0bZ0ZIRv3HBDhlHlz3ZrpRonTsAtt2QTlEw2oJlwzM6ebw0fOYWjHXdaVzdVKivcddcTN02lckRRvRJjvAd4YZu7XtTL8zbKZc7PzlIqlXIzdSpv/BLdmSd95CObXiOAZG2NJ33kI5x/8Ysziip/thspS+69t8+RaKOhn0alzjjtRd1ULje4445Vjhy5YNUuaYhd7Eu0HuN0s85sN1KWXnttnyPRRo5sqCNOe1G3HTyIo2LSkPNLdGecbtaZ7fZtSU+cyC4oObKhzriQXpLUbX6J7ky7ggPpoUMWHNhifa3UhSNHePjwYS4cOcLy3FyuyigPI0c21BFL30qSus0dxDuzseBAYXGRxsQEIydP0nAH8cdZXyul/Mgk2QghXAX8BHAcOBFj/PA2xxWA9wLPoVkC6c0xxi/1K05t5mZ1kqRu2voleqRcZvnYMat2tbH1S3SpVAILD2gPyGpkYwL4JnDbJY47ClwRY3x5COElwC3Adb0OTpIk9cfGL9GlUomGX6ClgZJJshFjrAEfCCGcuMShh4FPhhBeQ7Mm+rPaHRRCuAm4qfXcu95qfZi2jt+tvMUDxtSpvMWUt3jAmCRJ6ra9sGbjHcCngNcAf9bugBjjrcCtrZvpbrd/3+tbx/dD3uIBY+pU3mLKWzwwmDEdOHCgi9FIkrQzuapGFUI4GkL4+IamLwBfjDG+m+aGTF/OJjJJknamUK8zNjPDyHXXMTYzQ6FezzokaeCt/97tP3rU37ucyGqB+NOATwAHgIdDCD8QY3wdMAk8e8OhHwWuCyF8DlgDfrLvwUo7VK8XOH68QK22n8lJq3ZJw6hQrzM+Pf1olaV9QHF+nuW5ORc/Sz2y9fcOmr936ac/DVbuykxWazbuoTlSsbV9FpjdcHsVuKGPoUmXpV4vMD09Tq1WoFlADebni+6OLQ2Z0Wp10xcegGKtxmi1allO7UqhXm9W7VpYaO5PcvKkX6C32O73rnHiBNxySzZBaU+s2ZD2jGp1lFqtuKmtVitSrY66W7Y0RNwZW93U7op9etddFG67zZGyDbb7vUvuvbfPkWijXK3ZkPa6hYVC2/bFxfbtkgbTI9tccX7k6qv7HEn+bZxjX7jxRufYt9Huin1y9iyj1WpGEeXTdjvSp9de2+dItJEjG1IXjY4+0rb96qvbt0vSMHvcFfszZxg/c8a1LVs4UtaZ7XakT0+cyC4oObIhSVK3XbGy0r79wQf7HEm+XWxtix6z3RX7xsREnyPJt/Ud6S8cOcLDhw9z4cgRlufm4ODBrEMbao5sSF20stI+f3/wQfN6aZj45bAzXrHvTLsr9umhQ6xUKhlGlU8bd6RXPvgNSB2r1wvceGOBo0f3MzMzRr3uOoStJifbV5yamLASlTRMVioVVqemNrWtTk355XAL17Z0pt0V+9U77nCqmfYERzbUEUu6dqZSWWF+vripItXU1CqVSvspFZIG0/qXw9FqlauWl3lofJyVSsUvh9q1rVfsS6USLC1lGJHUGZMNdcSSrp0plxvMzS1z6lSJen2NiQk39ZOG1fqXw1KpxHm/FLbl2hZp8JlsqCOWdO1cudzg9OkGS0v3Zx2KJOWaa1ukweeaDZpThGZmxrjuuhHXImzDtQiSpG5zbYs0+IZ+ZOOxtQjrU4T2uRahjUplhc9/vsi5c49NpTpwwLUIkqTd27i2pbC4yEi5zPKxY65tkQbI0CcbrkXoXJIkF70tSdJObVz4XCqVaLi+RRooQz+NyrUInalWR/na1zbnpl/72gjVavuyhZIkSdLQj2y4FqEzJmXqtrvvhptvHmNhocDkpFW7JEkaREOfbLgvQmdMytRN9XqB668vcvbsEx5tc62UBk2hXme0WmVkeZkx99mQNKSGPtlY3xehWh1lefkqxscf8gprGyZl6qZqdZSzZzev+XGtlAZJoV5nfHqaYq0GwD6gOD/P8tycCYd2ZT15LSwsNEsGnzwJ2+zALuXJ0Ccb0Ew4ZmfPUyqVWFryi047blbXuXq9wPHjBWq1/U4P2obT8jToRqvVRxONdcVajdFqddMu0FIntiavAOldd1G47TaTV+WeyYY65mZ1l/ZYKeUC0Pzi7PSgx3NangZdYWGhffviYp8j0SBol7wmZ8+avGpPGPpqVFI3XayUsh5Tqaxw6FC6qc1peRok7oytbjJ51V5msiF1kdODOlMuN7jjjlWOHLnA4cMPc+TIBUd/NFDcGVvdZPKqvcxpVFIXOT2ocwcP4mJwDayNO2NftbzMQ1aj0mVYqVQozs9vXrNx6JDJq/YEkw2pi6zaJWnd+s7YpVKJ8+6KrcuwMXktLC7SmJhg5ORJGlaj0h5gsiF1kVW7JEm9sJ68riuVSmASqz3AZEPqMqt2SZK6zX02tFeZbEiSJOWY+2xoL7MalaRM3H03zMyMcfTofmZmxqjXrdglSe1cbJ8NKe8c2ZDUd/V6geuvL3L27BMebXPzQ0lqz302tJc5siGp76rVUc6eTTa1ufmhJLXnPhvay0w2JPWdmx9qGBTqdcZmZhi57jrGZmYo1OtZh5RL66/T/qNHKdx4o69TG+02iXSfDe3Wxt+5fnw2ZTKNKoRQAN4LPAcoAG+OMX5pm2MfAu5s3TwVY7y9P1FK6hU3P9Sg27qgdx9QnJ9neW7OBb0bPG7h85kzjJ854+u0hftsqFvaFRtY/2yiVOrJOXc9shFCuPPSR23rKHBFjPHlwM8Ct1zk2IUY4ytaPyYa0gCoVFY4dCjd1Obmh7qYy+xz+q7dgt5ireaC3i18nTq3vs/G/R//eHO/jYMHsw4pl/p91X6vyeJ37qIjGyGEv36Ru9tPIOzMYeCTIYTXAG8GnnWRY9dCCL8HLAFvjzGebRPnTcBNADHG5kY3uzAyMrLrx/ZK3mLKWzxgTJ3KU0ylEnzmMynveMcj3HtvwrXXppw4kXLw4DWZxpWn12hdHmPqlR72OX3ngt7O+Dqpm7a7ap9++tPuSdKSxe/cpaZR/XfgHiBpc9+1l3nudwCfAl4D/Nl2B8UYnwkQQvgR4JeAv9/mmFuBW1s306Vd7qhZKpXY7WN7JW8x5S0eMKZO5S2m7/iOErfcsjmerMPL22sElx/TgQMHuhhNz/Wyz+krF/R2xtdJ3bTdVfvGiRNwy8Um0QyPi/3O9WptxaWetxZj/Gvt7gghfOUyzvsFoBhjfHcI4fuAL7ee8yjw2hjj32vzmAT4y8s4pyQp33rV5/TdSqVCcX5+0xef1akpF/Ru4eukbtruqn1y7719jiS/LvY716u5BZdKNuJF7vsfl3HejwLXhRA+B6wBP9lqnwSevX5QCOEpwG8CDwP3AW+9jHNKkvKtV31O321c0HvV8jIPjY+zUqm46HmLrQufR8pllo8d83XSrmx31T69dk8NjPZUu2IDvf5sStI0vfRRe0t67ty5HT2gXi9QrY6yvHwV4+MPUams5GZjsbxN68hbPJC/mOr1AqdOlajV1picbOTm/ZS31ynD42KVAAAgAElEQVRv8cBgxtSaRtVuWpJ2Zsd9y7o8vq8gn3HlMSbIZ1zG9HiFep39R48y8rWvPdq29tSn8sjv/A5LOVuzkfVr1c5OYtpJ39Lx9KwQwpOA/RvbYox7fol/vV5genqcWq3YatnnTsbatcfeTwWaVZ3dGVvajUHtcyT11taL6AN4UX3P6aj0bQjhp4BzwP8CPg98Bbijh3H1TbU6uiHRaHInY+2W7yfp8g1ynyOpd0arVYpbRiCL585ROHEim4AEdL7Pxk8BU8BfxBivBd4I/H7PouojdzJWN/l+krpiYPscSb3jAvF86jTZSGKM54FHAGKMvw68smdR9ZE7GaubfD9JXTGwfY6k3nGBeD51umZjJYTwBOCuEMKvAv+zhzH1VaWywvx8cdPUF3cy1m75fpK6YmD7HEm9s11Z19RpVJnqdGTjeporzn+aZoLySuD1PYqpr8rlBnNzyxw5coGXv/wRjhy54GJe7dr6+2l6usHhww/7fpJ2Z2D7HEm9s17W9cKRIzx8+DAXjhxheW4ODh7MOrRcKdTrjM3MsP/oUcZmZijUe1t7o6ORjRjjl1t/vQ94U+/CyUa53GB29nyr5Nf5rMPRHlcuNzh9usHS0v1ZhyLtSYPe50jqnUa5zPnZ2azDyK1Cvc749PSm0Z/i/HwzKSuVenLOjpKNEMKHgMfVDosx/njXI5IkDbVB6XMK9Tqj1Sojy8uMuamfpBwYrVY3JRoAxVqN0WoV5uZ6cs5O12z80Zbbrwc+191QlHf1eoHjxwvUavtztVld3vg6SZdtz/c5W68e7uOxq4cmHJKyUvjqV9u3b0lAuqnTaVS/uvF2COHXgdt7EpFyyc3qOuPrJF2+fvQ5IYQC8F7gOTR/Wd8cY/xSt57/YlcPneIhKStXfP3r7dvvu+/xw8ndOucuH5cC7euLaSC5WV1nfJ2knuhFn3MUuCLG+HLgZ4Fbuvnk29X7LywudvM0krQj6VOesqP2buh0zcbX2Tx/dhT45Z5EpFxys7rO+DpJl69Pfc5h4JMhhNcAbwae1SaOm4CbAGKMlHaweLIwNQVnzjyufaRc3tHz9NLIyEhuYlmXx5ggn3EZU+fyGFdWMRW+8zthfr5te9KjmDpds/HCDX9PgftjjN/oejTKLTer64yvk9QV/epz3gF8CngN8Gdb74wx3grcuh7H0tJSx09cOHaM8TNnHlfvf/nYMRo7eJ5ealZgzEcs6/IYE+QzLmPqXB7jyiqmi302XbO21nFMBw4c6Picna7Z6N2qEe0JblbXGV8n6fL1qc/5AlCMMb47hPB9wJcv9YCdWK/3P1qtctXyMg9ZjUpSDmz8bCosLtKYmOj5Z9NFk40QwiO0KT/YaluJMV7Tk6iUO+ub1Z06VaJeX2NiwipL7fg6SbvX5z7no8B1IYTPAWvAT3bxuYHH6v2XSiXO5+yqqqTh1e+9SC41sjFKcxfXPwWeCzSAIs0KHv+xt6Epb9ysrjO+TtKu9a3PiTGuAjd08zklSY930WRjfY5sCOFCjPGB9fYQwhuAu4Bf7214kqRhYZ8jSYOn0wXi/y2E8C7g/cDDwHU0rzZJktRt9jmS1COFer25ZmNhgcbkZLZrNjY4RrMG+TzNjVDvwuFnSVJv2OdIUg8U6nXGp6c3VaMqzs+zPDcHPSrF22k1qvPAG3oSgSRJG9jnSFJvjFarmxINgGKtxmi1CnNzPTnnbncQJ4RQ6WYgkobL3XfDzMwYR4/uZ2ZmjHrdjQ+1PfscSbp8hYWF9u2Liz07Z6fTqNp5I1DtViCShke9XuD664ucPfuER9vm54vMzS1bJniDer3A8eMFarX9TE4OfRll+xxJukyNycn27RMTl5UUXMyl9tn4DeCrwMu23JUAT+tRTJIGXLU6ytmzyaa2Wq1ItTrK7Oz5jKLKl3q9wPT0OLVaAWiO+gx6QmafI0m9tVKpUJyff9wO4iuVCr3aPO9SScw54D7gWjYvzkuwBKGkXVpYaD9lanHRqVTrqtXRTTvRw1AkZPY50ja2VhDi5EkYHc06LO0xudtBPMb4/wCEEH4sxvi5jfeFEB7qWVSSBtrkZPsr8xMTg3nFfjeGMSGzz5Haa1dBKL3rLgq33dbTL4kaTP3eQbyjBeIxxme3af6hLsciaUhUKit8x3ekm9qe+tQ1KpWVjCLKn2FOyAalzyneeSff/qIXUXzKU5p/3nln1iFpj2pXQSg5e7ZZQUjKuV1Xo4ox/lk3A5E0XNJ06+20/YFDqlJZYWpqdVPb1NTq0CZke63PKd55J6XXvpbiPfeQ/NVfUbznnuZtEw7tQhYVhKRu6SjZCCG8OoRwRevvbw4h/NsQwnN6G1r/1OsFZmbGuO66EUtwSn1QrY5yzz2bF4ifO9dcj6CmcrnB3Nwy09MNDh9+mCNHLgz04vCNBqHPGTt2jGRtbVNbsrbG2LFjGUWkvexiFYSkvOu0ytWvxBg/GUL4TqAC/CLwPuClPYusTx6r+LK+EHPfwFd8kbI2jOsRdqNcbnD6dIOlpfuzDqXf9nyfU3jggR21SxfTroJQeugQKxW3n1H+dZpsfLP1548C/zLG+MEQws27PWkIoQC8F3gOzZqOb44xfmm3x12OIa34ImVqmNcjqCNd7XOy0Hjyk7miTWLRePKTM4hGe127CkIjJ0/SsBqV9oBOk43FEMLbgB/nsStLu17vARwFrogxvjyE8BLgFuC63R4XQrgJuAkgxkipVOo4kOXl9i/B8vJVO3qeXhkZGclFHOvyFg8YU6fyFNPJk3DXXemmvTYOHUo5eTLbGPP0Gq3LY0x90O0+p+/OnzpF6bWv3TSVKh0Z4fypUxlGpb1sawWhUqkES0sZRiR1ptNk403AzcDbY4y1EMLVwH+5jPMeBj4ZQngN8GbgWZdzXIzxVuDW1s10aQe/fOPjY8C+Nu0PsbSU/chGqVRiJ/+eXstbPGBMncpTTKOjcMcdJW6+eY3FxQITE83dsUdHG5n2nXl6jdZdbkwHDhzoYjR90+0+p+9WX/xilj72McaOHWNkZYW10VHOnzrF6otfnHVoktRXHSUbMcb/CfzYhtsPAjde5rnfAXwKeA1wsSojnR63K5XKCvPzxU1TqYa54ovULwcP4lTFS6jXCxw/XqBW28/kZDMhG4a1ZD3qc/pu9cUv5uuf/3wuk1hJ6peOh6VDCC8NIfxY6++FEMK3X8Z5vwB8Mcb4buCFwJdbz3s0hPDxSx3XTesVX44cucDLX/7IUFV8kZRf68Ur5uYKnDlzJbffvo/p6fGhqZbX5T5HkpSRjkY2Qgg/D/wdYBL4EJAAvwM8b5fn/ShwXQjhc8Aa8JOt9kng2R0c11XlcoPZ2fOtq09eaZWUvWEuXtGDPkeSlJFO12y8DvibwBcBYoxrIYTixR+yvRjjKnBDm/ZZYPZSx0nSoPvqV9uPYNRqQzGy0dU+R5KUnU6nUX0rxvjQ+o0QwihwVW9CkiR9/evtP57vu29PFWXaLfscSRoQnfZanwwhnAKeFEL4B8B/AH67d2FJ0nB7ylPSHbUPGPscSRoQnSYbPw/cBfw34EeA3wTe3qugJGnYTU2t7ah9wNjnSNKA6LT07SPAv2r9ABBCSLZ/hCTpcgxzWW77HEkaHBdNNkIIBeAfAE8H/kuM8Xda7dPAO9lcOUqS1CXrZblPnSpRr689uvHhIJflts+RpMFzqZGNW4FnAncC/zyE8GvANPAE4G09jk2SBKRDsUwDsM+RdJkK9Tqj1SqFhQUak5OsVCpQKmUd1lC7VLLxMuCvxxhXQwj/BFgE3hJj/EDvQ5Ok4bW+qV+z1G2z3O38fHHQNx0dqD5n/UvPyPIyY+PjrFQqNMrlrMOSBlahXmd8eppirfZoW3F+nvTTn4bR0QwjG26dLBAfCSHsAx4GasBtIYR9rTZJUg9cbFO/ATcQfc76l559t9/OFZ/7HPtuv53x6WkK9XrWoUkDa7Ra3ZRoABRrNQonTmQTkIBLj2w8A3iQ5u6t69Zvp6xfbpMkddXCQvuP18XFgf7YHZg+Z7svPaPVKudnZ7d5lLS9rdODOHnSq/VbFBYW2rYn997b50i00UWTjRjjUOweJUl5MznZfqrUxMTATqEaqD5nZEuical26WLaTQ9K77qLwm23OTVvg8bkZNv29Npr+xyJNhqYD3ZJGiSVygpTU6ub2oal9O0gSO67b0ft0sW0GylLzp5ltFrNKKJ8WqlUWJ2a2tS2OjVFw2lUmeponw1JUn8NY+nbQfLIt3873HPP49uf8pQMotFet930oMLiYp8jybdGuczy3FxzutniIo2JCVYqFa45eBCWlrIOb2iZbEhdVq8XOH68QK22n8lJvyBq98rlBqdPN1hauj/rULRDjac/Hb74xce3b7nqKnViu+lBjYmJPkeyhwxRzfC8M9mQumhIy5VK2mKlUqE4P79p6svq1FSz5r+0Q+3eT+mhQ76ftrD0bT65ZkPqoiEuVyppg/XpHBeOHOGRl7+cC0eOsDw352Je7crG99PDhw9z4cgRVu+4w/fTFpa+zSdHNqQuGtJypZLaaJTLnJ+dpVQqcd754rpM6++ndaVSyXUIW1j6Np8c2ZC6aBjLlUqSlAePbDNVKnUKVaZMNqQuslypJEnSY5xGJXWR5UolScrGFSvtL+wl27SrP0w2pC6zXKkkSf3nDuL55DQqSZIk7XnuIJ5PjmxIkiRpz3MH8Xwy2ZAkSdJA2FoiWNlzGpUkSZKknnBkQ5IkSQOhUK83p1EtLNCYnGSlUoFSKeuwhprJBnDnnUWOHRtjZWWE0dFv59Sp87z4xauXfqAkSVIfbP0SzcmT4GZ1mxTqdcanpynWao+2FefnST/9aV+rDA39NKo77yzy2teWuOeeIn/1Vwn33NO8feedxaxDkyRJevRL9L7bb+fKM2fYd/vtFF/1Kgr1etah5cpotbop0QAo1moUrEaVqaFPNo4dG2NtLdnUtraWcOzYWEYRSZIkPabdl+jk7FlGq9WMIsqnwsJC2/bk3nv7HIk26vs0qhDCq4GbgQLwsRjjP7vIsR8Gng+cB/48xvimbsfzwAOFHbVL6o6774abbx5jYaHA5KQ7rUvSdrb7El1YXOxzJPnmpn751NdkI4SwD6gCLwa+AfxxCOE3Y4wXGwf86RjjZ3sV05Of3OCBBx4/wPPkJ/ulR+qVer3A9dcXOXv2CY+2zc8XmZtbNuGQpC22+xLdmJjocyT5tlKpUJyf3zQKtDo1Reo0qkz1LNkIIbwf+O4tzceBPwW+Dfgg8DDwAmC7ZGMJ+KchhG8C748xfnSbc90E3AQQY6S0g6oDH/pQyg/+YLppKtXISMqHPpTu6Hl6ZWRkJBdxrMtbPGBMncpTTMePFzh7dvP0xVqtyKlTJU6fzi7ZyNNrtC6PMUnqr3ZfotNDh5qVlvQoN/XLpyRN076dLITwAiACX6KZeMwAvxdj/K1LPG4MmAeeG2P8xiVOk547d25HcW2uRrWWq2pUpVKJpRz9guQtHjCmTuUppqNH93PmzJWPaz98+GE+/vH7M4ioKU+v0brLjenAgQMAyaWO0yXtuG9Zl8f3FeQzrjzGBPmI69FqVK0v0SMnT7KUswpLeXid2sljXHs9pp30Lf1es/EVmovSbwDWgFcAtwCEEK4BPgG8Ncb4J1se16A5CvKtXgT14hev8vnPfz2X//HSIJqcbD96MTHhFCpJamfrztilUsmr9doT+ppsxBgfDCG8A/gssAq8L8a4fqnoicB3AY+WgQohnAKeSzPZeGuMMR/DDZIuS6Wywl13PXHTVKqpqVUqlZUMo5IkSd3W92pUrXUXj1t70Uo6Jra0HetXXJL6p1xu8P73r/JjP5bwwAMFnvzkBu95z3kXh0uSNGCGfp8NSf1Xrxd405uK3HNPkQceuIJ77inytreNUa9bclqSpF4q1OuMzcyw/+hRxmZmer45ZN9HNqRBV68XOH68QK223/0jtlGtjratRlWtjjI7ez6jqDSIQggvA34OuDLG+MpW25OAXwMO0Jyme2OM8Z7sopSk/ljfjX5jZbPi/DzLc3PQo8qHjmxIXVSvF5ieHmdursCZM1dy++37mJ4e94r9FgsL7V+PxUVfJ3Xd84BfYHPVlBngrhjjK4APAf8og7gkqe/a7UZfrNV6uhu9IxtSF1Wro9RqxU1tXrF/vNHRR9q2X311+3bpUrbZ2+kNMcbZEMLTt7QfBn4+hPAG4IeBp23znLvew2mjvO6Vkse48hgT5DMuY+pcHuPKKqaR5eW27VctL0OPYjLZkLrIK/ZSNmKMb9rhQ2aBjwHXA5/b5jlvBW5t3Ux3Wxo9r2XV8xhXHmOCfMZlTJ3LY1xZxTQ2Ps6+Nu0PjY8zsra20302OuI0KqmL3D+iMysr7T96HnzQjyT1xReAO2KM7wP+b+D3M45HkvpipVJhdWpqU9vq1FRPd6N3ZEPqokplhfn54qapVO4f8XgmZeqXEMIHgO8FDoYQ/gh4LXAK+EgI4VXAg8CPZxiiJPVNo1xmeW5u0270K5UKjXK5Z+c02ZC6qFxuMDe3zKlTJer1NSYmrEbVjpv6qV9ijG/c5q4f6msgkpQTW3ej7zWTDanLyuUGp083WFq6P+tQcqtcbnDHHavcfPMai4sFkzJJkgaUyYakTBw8iBW6JEkacK7GlCRJktQTjmxIysTdd8PNN4+xsFBwp3VJkgaUyYakvqvXC1x/fZGzZ5/waNv8fJG5uWUTDkmSBojTqCT1XbU6uqkSFTy207qk4VKo1xmbmWH/0aMUbryRQr2edUiSusiRDUl9507rnanXCxw/XqBW2+9UMw2kQr3O+PQ0xVqt2XDmDONnzrA8N9fTuv+S+sdkQ1LfuanfpdXrBaanx6nVCkAzCXOqmQbNaLX6WKLRUqzVGK1W+7oPgKTecRqVpL6rVFY4dCjd1OamfptVq6ObdqIHp5pp8BQWFtq3Ly72ORJJveLIhqS+c1O/S3OqmYZBY3KyffvERJ8jkdQrJhuSMuGmfhfnVDMNg5VKheL8/KapVKtTU6xUKhlGJambnEYlSTlUqawwNbW6qc2pZho0jXKZ5bk5Lhw5wsOHD9OYnnZxuDRgHNmQpBwqlxvMzS1z6lSJen3NqWYaWI1y+dHF4KVSicbSUsYRSeomkw2aVV+q1VGWl0cYHx+zQ9dlsVypuqVcbnD6dIOlpfuzDkWSpF0Z+mTjsfKS61Vf9lleUrtmuVJJkqTHDP2aDctLqpt8P0mSJD1m6JMNy0uqm3w/SZIkPWbop1FZXlLd5PtJknamUK8zWq1SWFigMDVF4dgxq1FJA2ToRzYqlRWe+tS1TW1Pfeqa5SW1K5YrlaTOFep1xqen2Xf77Vx55gyFuTnGp6cp1OtZhyapS4Y+2QBI0/Sit6VOrZcrnZ5ucPjwwxw5csHF4ZK0jdFqddOGfgDFWo3RajWjiCR129BPo6pWRzl3bvOC3nPnmgt63d1Yu2G5UknqTGFhoX374mKfI5HUK31PNkIILwN+DrgyxvjKSxz7auBmmjVEPxZj/GfdjscFvZIkZaMxOdm+fWKiz5FI6pUsplE9D/gFILnYQSGEfUAV+EHg+4E3hhC6vmLMBb2SJGVjpVJhdWpqU9vq1BQrlUpGEUnqtp6NbIQQ3g9895bmN8QYZ0MIT+/gKZ4F/CnwbcAHgYeBFwCPWzUWQrgJuAkgxkipVOo4zpMn4a67Us6efSz3OXQo5eTJkR09T6+MjOQjjnV5iweMqVN5iylv8YAxSf3WKJdZnptrVqNaXGSkXGbZalRST22sANeYnGSlUunp71zPko0Y45u68DTPB04Bx4GZi5zrVuDW1s10aWmp4xOMjsJttxWoVkdZXr6K8fGHqFRWGB1tsIOn6ZlSqcRO/j29lrd4wJg6lbeY8hYPDGZMBw4c6GI0Uvc1ymXOz84Czfd7I2e/g9IgWa8At7EwQ3F+nuW5OejRha3cLBAPIVwDfAJ4a4zxT4Cv0JzmdQOwBrwCuKUX5y6XG8zOnm916i4KlyRJ0uC5aAW4ubmenLPvazZCCB8A/h3w3SGEPwohPKN11xOB7wLGAGKMDwLvAD4L/GfgfTHGc/2OV5IkSRoEWVSA6/vIRozxjdu0nwMmtrR9FPhoP+KSJEmSBtnFKsD1KilwUz9JkiRpCGRRAS43azYkSZIk9c7WCnCNiYm9W41KkiRJUr5srADXD06jkqScqtcL3HhjgaNH9zMzM0a9Xsg6JEmSdsSRDUnKoXq9wPT0OLVaAWgmGfPzRebmlimXG9kGJ0lShxzZkKQcqlZHqdWKm9pqtSLV6mhGEUmSBkGhXmdsZob9R48yNjNDoV7v6fkc2ZCkHFpYaD9lanHRqVSSpN3JYgdxRzYkKYcmJ9tPlZqYcAqVJGl3LrqDeI+YbEhSDlUqK0xNrW5qm5papVJZySgiSdJeNxQ7iEuSLq1cbjA3t8ypUyXq9TUmJhpUKisuDpck7VoWO4ibbEhSTpXLDU6fbrC0dH/WoUiSBsBKpUJxfn7TVKr1HcSv6dE5TTYkSZKkIeAO4pKkR9XrBY4fL1Cr7Wdy0mlUkqTL1+8dxE02JCmH3NRPkjQIrEYlSTnkpn6SpEFgsiFJOeSmfpKkQWCyIUk55KZ+kqRBYLIhSTnkpn6SpEHgAnFJyiE39ZMkDQKTDUnKKTf1kyTtdU6jkiRJktQTJhuSJEmSesJkQ5IkSVJPmGxIkiRJ6gmTDUmSJEk9YbIhSZIkqSdMNiRJkiT1hMmGJEmSpJ5wUz9J0sAKIbwRuB4YB94fY/wXIYQnAb8GHAAawI0xxnsyDFOSBlbfk40QwsuAnwOujDG+8hLHfhh4PnAe+PMY45t6H6EkaYD8LvBB4GrgL4B/AcwAd8UYfzSEcAPwj4CbsgtRkgZXFiMbzwN+AfilDo//6RjjZy92QAjhJlodRYyRUqm0q8BGRkZ2/dheyVtMeYsHjKlTeYspb/GAMe1lIYT3A9+9pfkNMca7Wvc/Bai32g8DPx9CeAPww8DTtnnOge1bIJ9x5TEmyGdcxtS5PMY1TDH1LNm4yAf/bAjh6R0+zRLwT0MI36Q5/P3RdgfFGG8Fbm3dTJeWlnYTMqVSid0+tlfyFFO9XuDUqRK1WsrkZINKZYVyuZF1WLl6jdYZ06XlLR4YzJgOHDjQxWjy62Ij3yGECeBfAT++oXkW+BjNKVaf2+Y5B7ZvgXzGlceYIJ9xGVPn8hjXXo9pJ31Lz5KNbkx5ijG+HSCEMAbMhxD+XYzxG5cdnHasXi8wPT1OrVYACgDMzxeZm1vORcIhSe2EEA7RTDTeEmP801bzF4CHYozvCyH8feD3MwtQkgZcbhaIhxCuAT4BvDXG+Cdb7m4ADwPf6ntgAqBaHaVWK25qq9WKVKujzM6ezygqSbqk36LZ1/1qCAHg3cAp4CMhhFcBD7J5xEOS1EVZLBD/APC9wMEQwh8Br40x/gXwROC7gLENx54Cnksz2XhrjHG13/GqaWGh0LZ9cbF9uyTlQYzxudvc9UN9DUSScqJQrzNarVJYWKAxOclKpUKjXO7Z+fqebMQY37hN+zlgYkvbsb4EpUuanGw/VWpiwilUUq/U6wWOHy9Qq+3P1TopSdLeVKjXGZ+eplirPdpWnJ9neW4OerRgPTfTqJRvlcoK8/PFTVOppqZWqVRWMoxKGlyuk5Ikddtotbop0QAo1mqMVqswN9eTc7qDuDpSLjeYm1tmerrB4cMPc+TIBb/0SD10sXVSkiTtRmFhoX374mLPzunIhjpWLjc4fbrB0tL9WYciDTzXSUmSuq0xOdm+fWKiZ0mBIxuSlEOuk5IkddtKpcLq1NSmttWpKVYqlZ6d02RDknKoUllhampzAT7XSUmSLkejXGZ5bo4LR47w8OHDXDhyhOW5ucGqRiVJurT1dVKnTpWo19eYmLAalSTp8jXKZc7PzvbtfCYbkpRTrpOSJO11TqOSJEmS1BMmG5IkSZJ6wmRDkiRJUk+YbEiSJEnqCZMNScqper3AjTcWOHp0PzMzY9TrbugnSdpbrEaljtXrBY4fL1Cr7Wdy0jKcUi/V6wWmp8ep1QpAM8mYny8yN7fs750kac8w2VBH/OIj9Ve1OkqtVtzUVqsVqVZHmZ09n1FUkqS9rlCvM1qtUlhYoDE5yUql4qZ+yp5ffKT+WlhoP2VqcdGpVJKk3SnU64xPT1Os1R5tK87Pszw3B6VST87pmg11xC8+Un9NTrYfMZyYcCRRkrQ7o9XqpkQDoFirMVqt9uycJhvqiF98pP6qVFaYmlrd1DY1tUqlspJRRJKkva6wsNC+fXGxZ+c02VBH/OIj9Ve53GBubpnp6QaHDz/MkSMXXCMlSbosjcnJ9u0TEz07p2s21JH1Lz6nTpWo19eYmLAaldRr5XKD06cbLC3dn3UokqQBsFKpUJyf3zSVanVqipVKhWt6dE6TDXXMLz6SJEl7V6NcZnlurlmNanGRxsSE1agkSZIkdUejXOb87GzfzueaDUmSJEk9YbIhSZIkqSdMNiRJkiT1hMmGJEmSpJ4w2ZAkSZLUEyYbkiRJknrCZEOSJElST/R9n40QwhuB64Fx4P0xxn9xkWNfDdwMFICPxRj/WX+ilCRJknS5shjZ+F3gB4CXACe2OyiEsA+oAj8IfD/wxhBC77Y3lCRJktRVSZqmPXniEML7ge/e0vyGGONdrfufQXO04oXbPP4FwM8BbwPeAzwDeFeM8bfaHHsTcBNAjHHrOSVp2CVZBzAAetNZStLe1VHf0rORjRjjm2KML9zys55oTAD/CvjxSzzN84FTQAX4zxc5163r56D5D9/VTwjhC5fz+F785C2mvMVjTHs3przFM+Ax6fJl/X+Y1/fWwMeU17iMaW/HNSAxdSSLNRuHaCYab4kx/umG9muATwBvjTH+CfAVmsnQDcAa8Argln7HK0mSJI4KJfUAAAhpSURBVGl3sliz8VvABPCrIYTPhhB+oNX+ROC7gDGAGOODwDuAz9Ic1XhfjPFc/8OVJEmStBt9H9mIMT53m/ZzNJOQjW0fBT7aj7habu3juTqVt5jyFg8YU6fyFlPe4gFjUm/k9f8wj3HlMSbIZ1zG1Lk8xjU0MfVsgbgkSZKk4eamfpIkSZJ6wmRDkiRJUk+YbEiSJEnqib4vEM+j1m7ltwJPBfYBr4sx/nm2UUEIYQb458AzY4xfzTgcQgjfC7wLuBK4Lcb4wYzjuQr418AB4CrgH8cY/21GcfwEcBw4EWP8cOu1+mWgAJyJMb49BzG9Gng7MAp8Ksb4jqxjarVfDfx34L/GGF+fZTwhhBHg3cD3A48AfzvG2Mg4ph+iWZnvW8BfAq+NMX6zXzFp9/Lat0D++hewj+kgFvuZXcbVas+kr9kupmHqbxzZAGKMF4B3xhhfCXwQeEvGIRFCOAj8EPAHWccCEEIYpbnB4utijD+QdSfQ8iLg4RjjS4A3Aa/PKI4J4JvAbRvabgVuiDF+P/D8EMLhHMT0x8D/AXwv8LrW3jZZxwTwT4Df6HMs0D6efwgsxhhfHmN8ZT8/+C8S09uAH40xvpTmJkrP7nNM2qU89i2Qv/4F7GM6YD9zeXFBdn0NDHl/Y7LREmP8i9ZfDwBns4wlhJAAvwL8FM1sNw9eCjwBOB1C+C8hhNdkHRBwJ3B1COEDND9EfjGLIGKMtRjjB4BVeHSDylXgwRD+//buPUaPqg7j+LcWa6DtKoiCSmyIKxivXL0QgSaoMQ2UxtgnxEuzSmmK0YQmJooQU6tNEaoYpGqLVrmI9DGmXgioESzWEkGNSZUIhYR4ISkBQ1IqtIDWP86pLCu7fXdfZ2e6+3z+2vfMOzO/3bw7T86ZM+/Rt4A5wKlt1lTb/mr7acqI0x5gd9s1STqTMor508msZbR6gHOBd0i6TdL1kmZ3oKa1wPWS1gMPAX+YzJqiP13KFuhsvkAyZkzJmf7qajNrRquJaZQ36WwMI2kRcCKwruVSlgG3297Rch3DDQC/s70AeB9lBKptA8Bs4G7KP8ukXmgP4JXADcBV/O/oSmvqxexG4GM1ENqs5TBgJeWWe1cMACtsnwXsBJa0XA/APGAHcC/l+tTGSGH0oUPZAt3MF0jGTERyprd6upg1MI3yJp2NStL5wGJgse1nWi7nHGCRpC3ACcBNko5vtyTup3wIoczl+2eLtey3ENhuewMlQC9ouR4AbD9GGc1ZSrmlvIAOTFeQdCSwGfii7V+0XQ9wGvBS4EeUuePvlbSm3ZKe8znfS7nF3LZPA5+wfSXwR+CMluuJcehYtkA38wWSMeOSnBmXLmYNTKO8yaJ+gKRTgLuAbZTbyk/Zfk+7VRU1EIa68ACfpK9Q5mE+DayyfVvL9byCMprzAsqXHVxp+wct1HEM8EPKKNNeYCvwHWAN5W/1c9urOlDT4cAbgb/Ut22wfWObNdleUrfNp3zOh9qsB/gy8DXK7f+/AUsnc2RulJp+Sxnx2g08CnzE9qRPTYjx63K2QLfyBZIxB6glOdNHXW1mzWg1MY3yJp2NiIiIiIhoRKZRRUREREREI9LZiIiIiIiIRqSzERERERERjUhnIyIiIiIiGpHORkRERERENOKQtguIGK/6dY3HAU/Upm/avmzY9p8AD9heMWK/y4FbbG8Z0X4kcF1dTGoi9XwcOKWJr9KTdDKwCTgC2Gi7a4sSRURMCcmWiGaksxEHq2W2bx5l24PA35+n/a3Ar56nfU7d1jm2fw8MSlpJqTMiIpqTbIn4P0tnI6YMSUcDvwbmUhan2d++mLLw0auADZKeAHbZPknSTcDbgcMlPVB3+bbt1ZIOBb5KWX10JrDW9jX1mO8ErgBeRhkZ+vEBansDZXXXo2zvkTSDElwfsH2npE3AycAM4B7gg7YfP8Ax9wFzbe+WdDbwSdvzJc0EVlNWCp5FGVn7fN3neGA9cDRwKGWl0DFrj4iYzpItyZboT57ZiCnD9k7bg8CnRrR/v7bfRRm1GrR9Ut12HjAfeKy2D9peXXf9DHCf7ddTQuMSScdKegnl9vPSetzP9VDbPcC9wMLadAbwpO076+vl9dyvAR6nrOA5UUPAbMoqrm8GFko6s267GLjZ9utszwNu6eM8ERFTXrLlv4ZItsQE5M5GHKy+Lmlt/fli25sbOMcC4AhJ59fXs4DXUkai/lwv8gD/6vF4G4EPAwY+VF/v91FJS+o5Xg7c32fdbwPeXV/PAQaBO4BbgS9JGgC+a/u+Ps4TETHVJFvGrjvZEuOWzkYcrC4cY17taPaNs30mMGT7juGNkt4P7B3W1Ov/0feAyyUdA5wLvKke7zRgBXCi7UckXTbGMXoxE7jE9rUjN9jeJOlu4DzgZ5KusL2uz/NFREwVyZbRJVtiQjKNKqaTR4C3ANR5rfv9AxiQ9OoR27YCF0l6UW0/rLZvB06V9GJJxwEX9XJy27so82/XA9tsP1w3HVVre1TSIM/eDu/l95knaRawfFj7VmB5HWFC0gslHVJ/nmv7QdtrgC/w7AhVRERMTLIl2RJjyJ2NmDIkvQv4BuUhvtn1obyVtm+ob1kDXCdpKbCLGg71IbhLgW2SngKuBVYBnwWuBnZIehLYA5xge4ek9cCfKA/iXUWZv9qLjcDtPPeifyuwrB7rIeA3w36nIeBSyoOCMyQtAi6w/UtgJbAZeBjYApxed7saOBbYLmkP8AxwVn3fOkmnA/8GdgIX9lh3RMS0lGxJtkR/ZuzbN9pdvoiIiIiIiInLNKqIiIiIiGhEOhsREREREdGIdDYiIiIiIqIR6WxEREREREQj0tmIiIiIiIhGpLMRERERERGNSGcjIiIiIiIa8R/J4OGRIMrsPwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 估计值\n",
    "homo_x = homo_fitted_result.fittedvalues \n",
    "hetero_x = hetero_fitted_result.fittedvalues \n",
    "\n",
    "# 残差值\n",
    "homo_e = homo_fitted_result.resid\n",
    "hetero_e = hetero_fitted_result.resid\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 6))\n",
    "fig.tight_layout(pad=5)\n",
    "\n",
    "axes[0].scatter(homo_x, homo_e, color='b')\n",
    "axes[0].set_xlabel(\"Fitted values\")\n",
    "axes[0].set_ylabel(\"Residual\")\n",
    "axes[0].set_title(\"Residuals vs. Fits Plot\")\n",
    "\n",
    "axes[1].scatter(hetero_x, hetero_e, color='r')\n",
    "axes[1].set_xlabel(\"Fitted values\")\n",
    "axes[1].set_ylabel(\"Residual\")\n",
    "axes[1].set_title(\"Residuals vs. Fits Plot\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "表1\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>均值(homo)</th>\n",
       "      <th>方差(homo)</th>\n",
       "      <th>均值(hetero)</th>\n",
       "      <th>方差(hetero)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>x=1</th>\n",
       "      <td>0.234120</td>\n",
       "      <td>0.795501</td>\n",
       "      <td>1.481288</td>\n",
       "      <td>28.638050</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=3</th>\n",
       "      <td>0.562194</td>\n",
       "      <td>1.467287</td>\n",
       "      <td>4.780266</td>\n",
       "      <td>93.906343</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=5</th>\n",
       "      <td>0.381242</td>\n",
       "      <td>0.531841</td>\n",
       "      <td>4.371954</td>\n",
       "      <td>53.184129</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=7</th>\n",
       "      <td>0.170099</td>\n",
       "      <td>0.637128</td>\n",
       "      <td>2.948230</td>\n",
       "      <td>91.746470</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=9</th>\n",
       "      <td>-0.557512</td>\n",
       "      <td>0.810359</td>\n",
       "      <td>-6.479940</td>\n",
       "      <td>158.830340</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     均值(homo)  方差(homo)  均值(hetero)  方差(hetero)\n",
       "x=1  0.234120  0.795501    1.481288   28.638050\n",
       "x=3  0.562194  1.467287    4.780266   93.906343\n",
       "x=5  0.381242  0.531841    4.371954   53.184129\n",
       "x=7  0.170099  0.637128    2.948230   91.746470\n",
       "x=9 -0.557512  0.810359   -6.479940  158.830340"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "homo_nd = homo_e.reshape(5, 10)\n",
    "homo_df = pd.DataFrame(homo_nd) \n",
    "homo_mean = homo_df.apply(np.mean, axis=1)\n",
    "homo_var = homo_df.apply(np.var, axis=1)\n",
    "\n",
    "hetero_nd = hetero_e.reshape(5, 10)\n",
    "hetero_df = pd.DataFrame(hetero_nd) \n",
    "hetero_mean = hetero_df.apply(np.mean, axis=1)\n",
    "hetero_var = hetero_df.apply(np.var, axis=1)\n",
    "\n",
    "print(\"表1\")\n",
    "pd.DataFrame(\n",
    "    data=np.vstack((homo_mean, homo_var, hetero_mean, hetero_var)).T,\n",
    "    index=[\"x=1\", \"x=3\", \"x=5\", \"x=7\", \"x=9\"],\n",
    "    columns=[\"均值(homo)\", \"方差(homo)\", \"均值(hetero)\", \"方差(hetero)\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如果拟合模型接近总体模型, 残差的均值应该为0, 方差相同.\n",
    "\n",
    "因为模拟数据中, 将误差故意随着X的变化, 拟合后均值接近0, 但是方差明显递增, 那么这种情况是不直接使用OLS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "表2\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x=1</th>\n",
       "      <th>x=3</th>\n",
       "      <th>x=5</th>\n",
       "      <th>x=7</th>\n",
       "      <th>x=9</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>x=1</th>\n",
       "      <td>0.795501</td>\n",
       "      <td>0.030915</td>\n",
       "      <td>0.226428</td>\n",
       "      <td>-0.447775</td>\n",
       "      <td>0.059800</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=3</th>\n",
       "      <td>0.030915</td>\n",
       "      <td>1.467287</td>\n",
       "      <td>0.654936</td>\n",
       "      <td>0.107649</td>\n",
       "      <td>0.399611</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=5</th>\n",
       "      <td>0.226428</td>\n",
       "      <td>0.654936</td>\n",
       "      <td>0.531841</td>\n",
       "      <td>-0.119510</td>\n",
       "      <td>0.099588</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=7</th>\n",
       "      <td>-0.447775</td>\n",
       "      <td>0.107649</td>\n",
       "      <td>-0.119510</td>\n",
       "      <td>0.637128</td>\n",
       "      <td>0.068148</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=9</th>\n",
       "      <td>0.059800</td>\n",
       "      <td>0.399611</td>\n",
       "      <td>0.099588</td>\n",
       "      <td>0.068148</td>\n",
       "      <td>0.810359</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          x=1       x=3       x=5       x=7       x=9\n",
       "x=1  0.795501  0.030915  0.226428 -0.447775  0.059800\n",
       "x=3  0.030915  1.467287  0.654936  0.107649  0.399611\n",
       "x=5  0.226428  0.654936  0.531841 -0.119510  0.099588\n",
       "x=7 -0.447775  0.107649 -0.119510  0.637128  0.068148\n",
       "x=9  0.059800  0.399611  0.099588  0.068148  0.810359"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 同方差情景: 协方差\n",
    "print(\"表2\")\n",
    "pd.DataFrame(np.cov(homo_nd, rowvar=True, bias=True),\n",
    "            index=[\"x=1\", \"x=3\", \"x=5\", \"x=7\", \"x=9\"],\n",
    "            columns=[\"x=1\", \"x=3\", \"x=5\", \"x=7\", \"x=9\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "表3\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x=1</th>\n",
       "      <th>x=3</th>\n",
       "      <th>x=5</th>\n",
       "      <th>x=7</th>\n",
       "      <th>x=9</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>x=1</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.028615</td>\n",
       "      <td>0.348112</td>\n",
       "      <td>-0.628964</td>\n",
       "      <td>0.074481</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=3</th>\n",
       "      <td>0.028615</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.741396</td>\n",
       "      <td>0.111336</td>\n",
       "      <td>0.366472</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=5</th>\n",
       "      <td>0.348112</td>\n",
       "      <td>0.741396</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.205306</td>\n",
       "      <td>0.151697</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=7</th>\n",
       "      <td>-0.628964</td>\n",
       "      <td>0.111336</td>\n",
       "      <td>-0.205306</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.094842</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=9</th>\n",
       "      <td>0.074481</td>\n",
       "      <td>0.366472</td>\n",
       "      <td>0.151697</td>\n",
       "      <td>0.094842</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          x=1       x=3       x=5       x=7       x=9\n",
       "x=1  1.000000  0.028615  0.348112 -0.628964  0.074481\n",
       "x=3  0.028615  1.000000  0.741396  0.111336  0.366472\n",
       "x=5  0.348112  0.741396  1.000000 -0.205306  0.151697\n",
       "x=7 -0.628964  0.111336 -0.205306  1.000000  0.094842\n",
       "x=9  0.074481  0.366472  0.151697  0.094842  1.000000"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 同方差情景: 相关系数\n",
    "print(\"表3\")\n",
    "pd.DataFrame(np.corrcoef(homo_nd, rowvar=True, bias=True),\n",
    "            index=[\"x=1\", \"x=3\", \"x=5\", \"x=7\", \"x=9\"],\n",
    "            columns=[\"x=1\", \"x=3\", \"x=5\", \"x=7\", \"x=9\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "表4\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x=1</th>\n",
       "      <th>x=3</th>\n",
       "      <th>x=5</th>\n",
       "      <th>x=7</th>\n",
       "      <th>x=9</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>x=1</th>\n",
       "      <td>28.638050</td>\n",
       "      <td>1.483913</td>\n",
       "      <td>13.585688</td>\n",
       "      <td>-32.239803</td>\n",
       "      <td>5.023214</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=3</th>\n",
       "      <td>1.483913</td>\n",
       "      <td>93.906343</td>\n",
       "      <td>52.394847</td>\n",
       "      <td>10.334264</td>\n",
       "      <td>44.756439</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=5</th>\n",
       "      <td>13.585688</td>\n",
       "      <td>52.394847</td>\n",
       "      <td>53.184129</td>\n",
       "      <td>-14.341250</td>\n",
       "      <td>13.942340</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=7</th>\n",
       "      <td>-32.239803</td>\n",
       "      <td>10.334264</td>\n",
       "      <td>-14.341250</td>\n",
       "      <td>91.746470</td>\n",
       "      <td>11.448874</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=9</th>\n",
       "      <td>5.023214</td>\n",
       "      <td>44.756439</td>\n",
       "      <td>13.942340</td>\n",
       "      <td>11.448874</td>\n",
       "      <td>158.830340</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           x=1        x=3        x=5        x=7         x=9\n",
       "x=1  28.638050   1.483913  13.585688 -32.239803    5.023214\n",
       "x=3   1.483913  93.906343  52.394847  10.334264   44.756439\n",
       "x=5  13.585688  52.394847  53.184129 -14.341250   13.942340\n",
       "x=7 -32.239803  10.334264 -14.341250  91.746470   11.448874\n",
       "x=9   5.023214  44.756439  13.942340  11.448874  158.830340"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 异方差情景: 协方差\n",
    "print(\"表4\")\n",
    "pd.DataFrame(np.cov(hetero_nd, rowvar=True, bias=True),\n",
    "            index=[\"x=1\", \"x=3\", \"x=5\", \"x=7\", \"x=9\"],\n",
    "            columns=[\"x=1\", \"x=3\", \"x=5\", \"x=7\", \"x=9\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "表5\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x=1</th>\n",
       "      <th>x=3</th>\n",
       "      <th>x=5</th>\n",
       "      <th>x=7</th>\n",
       "      <th>x=9</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>x=1</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.028615</td>\n",
       "      <td>0.348112</td>\n",
       "      <td>-0.628964</td>\n",
       "      <td>0.074481</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=3</th>\n",
       "      <td>0.028615</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.741396</td>\n",
       "      <td>0.111336</td>\n",
       "      <td>0.366472</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=5</th>\n",
       "      <td>0.348112</td>\n",
       "      <td>0.741396</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.205306</td>\n",
       "      <td>0.151697</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=7</th>\n",
       "      <td>-0.628964</td>\n",
       "      <td>0.111336</td>\n",
       "      <td>-0.205306</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.094842</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=9</th>\n",
       "      <td>0.074481</td>\n",
       "      <td>0.366472</td>\n",
       "      <td>0.151697</td>\n",
       "      <td>0.094842</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          x=1       x=3       x=5       x=7       x=9\n",
       "x=1  1.000000  0.028615  0.348112 -0.628964  0.074481\n",
       "x=3  0.028615  1.000000  0.741396  0.111336  0.366472\n",
       "x=5  0.348112  0.741396  1.000000 -0.205306  0.151697\n",
       "x=7 -0.628964  0.111336 -0.205306  1.000000  0.094842\n",
       "x=9  0.074481  0.366472  0.151697  0.094842  1.000000"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 异方差情景: 相关系数\n",
    "print(\"表5\")\n",
    "pd.DataFrame(np.corrcoef(hetero_nd, rowvar=True, bias=True),\n",
    "            index=[\"x=1\", \"x=3\", \"x=5\", \"x=7\", \"x=9\"],\n",
    "            columns=[\"x=1\", \"x=3\", \"x=5\", \"x=7\", \"x=9\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由上面的表可以看到, 对角线上的数据就是X在个点的样本残差的方差.\n",
    "\n",
    "除主对角线, 其他值基本都在0附近.\n",
    "\n",
    "对比同方差和异方差两种场景, 得到如下结论:\n",
    "\n",
    "1. 表2, 主对角线数值基本相等; 表4, 主对角线递增; 两个表的主对角线均和表1里面的数值有对应关系\n",
    "\n",
    "2. 表3,表5是相关系数表, 可以看出5组样本采集的数据是非线性相关.\n",
    "\n",
    "3. 表2和表4也是特殊的残差的协方差矩, 特殊在我X取5个点, 设定X后, 再抽出10个样本, 真实场景X的值不是如此设定的."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 显示统计信息"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.992\n",
      "Model:                            OLS   Adj. R-squared:                  0.992\n",
      "Method:                 Least Squares   F-statistic:                     6400.\n",
      "Date:                Sat, 29 Dec 2018   Prob (F-statistic):           1.36e-53\n",
      "Time:                        18:16:58   Log-Likelihood:                -71.433\n",
      "No. Observations:                  50   AIC:                             144.9\n",
      "Df Residuals:                      49   BIC:                             146.8\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1             2.0088      0.025     79.998      0.000       1.958       2.059\n",
      "==============================================================================\n",
      "Omnibus:                        3.737   Durbin-Watson:                   1.772\n",
      "Prob(Omnibus):                  0.154   Jarque-Bera (JB):                1.920\n",
      "Skew:                          -0.171   Prob(JB):                        0.383\n",
      "Kurtosis:                       2.103   Cond. No.                         1.00\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "print(homo_fitted_result.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.553\n",
      "Model:                            OLS   Adj. R-squared:                  0.544\n",
      "Method:                 Least Squares   F-statistic:                     60.60\n",
      "Date:                Sat, 29 Dec 2018   Prob (F-statistic):           4.08e-10\n",
      "Time:                        18:16:58   Log-Likelihood:                -187.11\n",
      "No. Observations:                  50   AIC:                             376.2\n",
      "Df Residuals:                      49   BIC:                             378.1\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1             1.9764      0.254      7.785      0.000       1.466       2.487\n",
      "==============================================================================\n",
      "Omnibus:                        1.333   Durbin-Watson:                   1.726\n",
      "Prob(Omnibus):                  0.514   Jarque-Bera (JB):                1.243\n",
      "Skew:                          -0.366   Prob(JB):                        0.537\n",
      "Kurtosis:                       2.756   Cond. No.                         1.00\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "print(hetero_fitted_result.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ---------------\n",
    " \n",
    " 异方差的$R^2$远小于同方差的$R^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "homo_fitted_result = sm.WLS(endog=y1, exog=x).fit()\n",
    "hetero_fitted_result = sm.WLS(endog=y2, exog=x).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Residuals vs. Fits Plot')"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxsAAAF0CAYAAABL4yRvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3X98pFd92PvP49Fgs1hE1k6QvJARu9AQWqDQkEA2/Myr1/eGS2g23R7U2MUkEBNAuZvLkmmccMtCkm47r5h2iVKKA4VtcBGHpCYNNYU2KaRN1jRB1CEpSWnWzNSsJSMrG8ssNtL4uX/MyJbk0e5IOzPPo5nP+/XSa3fOPDPPd2dHc+b7nHO+J0nTFEmSJEnqtiuyDkCSJEnSYDLZkCRJktQTJhuSJEmSesJkQ5IkSVJPmGxIkiRJ6gmTDUmSJEk9YbKhnkiS5ElJktydJMlsB8d+NkmSE1067wuTJEk3/Dy9G8+7lyRJciRJkvNJkrx0m/tHkyQpZhDX7Ib/l8/u8jmuSJJkrMuhSdoj7FuyY9+i3TLZGEKtD+A0SZJHkiS5N0mSDyZJMtrl06wBdWCxy897UWma/lGapglwsJ/n7Zctnd36z3dtOeyvgP8NPNjm8U8D/hJ4ag9i++yGmO5LkuQ3kiR55vr9aZrOtP5v3nUZp/kU8NOXHaykrrNv2bvsW+xbeslkY3i9CygArwC+D3hnN588TdOH0zR9eZqmv9DN5xUAr0zTNNnw82cb70zT9HfTNH1umqZfbPPYEZr/773yrtaH/ncDDeA/JEnyhC4+/5VdfC5J3WffsnfZt6gnTDaGWNr05zQz+kevYCRJUkyS5J8mSXJ/kiRLSZLcliRJacP91yRJ8q9bw6mrSZJ8LUmS79xw/yc3XIX48NbzJkny9CRJPpEkyTeSJFkBvnfDfa9IkiTdcvxXkyR5fevvT0yS5D8kSfL1JEkeTpLkfyRJ8iM7+XcnSfKPkiQ5lyTJWuvfeEOHjwutc16zoe1g6yre97Rul5Mk+a0kSR5MkuRbSZLUkiS5eifx7VaSJKUtV6VeseG+p7de17tbTXdvc9ybWlMUVlv/vz+7m1jSNP3fwAzwDOBvdhj/05Mk+e0kSS4kSXJPkiS/miTJk1v3fbYV/8uBd17s/SUpW/Yt9i32LdrIZGOIJUlSSJLk+cCPAP9xw10ngZcAfwt4HnAN8Gsb7n8nUAaeDTwJeBlwbv3ONE1f3boCcbrNOa8Cfgf4MvA0mh8Yf76DsFPgw624xlp//zedfugmSfLDwE8Bfxt4IvCcVjyd+ATwADC9oe11wF1pmv5h6/Z7gQs0X59R4AeBhzp8/suSpunS+hWpNvd9dcsUgIMbrl59FiBJkhcAs8CP0vx//WvAv7mMkEZaf379Uge23he/S/O98B3A/0nzPfjeVvyvaMX/OVpXuFo/r7+M+CT1gH2LfYt9izYy2Rhe76Q59/UPgVNpmp6C5iIp4M3Az6ZpWkvT9Fzr2B9KHhuyfJDmvMxnA2tpmv5FmqaPm8O5jR+h+UH882ma/mWapvfRnAfakTRNH0rTdC5N03vTNP0m8EGaw5+HOnyKB4Grge8Biq3nubfDc38L+NfAjQBJkiQ0O4Rbtzz/M4BntIb7/0eapmsdxtap/7zh6kunr3snLrT+fDGwL03Tr6dpWt/pkyRNT6fZudyRpulXO3jYq2m+L/5hmqb3p2n6p8DNwHTrPSlpb7BvsW/Zyr5lyPlCD693ASXgszTn1a4rAfuA31v/0AH+G825mJOtY34RmAM+AnwtSZJ/nCTJvg7P+3Tgq2maPrKboJMkuTJJkl9pDYWuAUutuzp6L6dp+p+AnwSOA19PkuTXkySZ2kEIHwBelCTJs4CX0nxNbttw/wzweeDTSZKcTZLkZ5Ik6fY81o3zars2jN6a9vDDQAAWW9MRnrvDp3knzfm0d9JcwPmjHT7uEM33RWND21manf3+HcYgKTv2LfYtm9i3yGRjiKVpej/wFuBokiQvaTXfD6wC37dloViyfiWidQXo/6N5BepHgb9P51UgltjwC54kyZVsrl6x1movbPjzqg33vw34v2guPizSHIbfkTRNT6dp+jyaC80mgN/YwWO/DPw+zatOrwM+lqbpAxvuP5+m6U8BT6FZ2eLnaHZAmyRJcqAHHUUn1j9wR9rdmabpv0/T9PuAZ9K8kvafWlfZOvWuNE2vSNN0Mk3Tt6Rp2umVxXuBqS1Xmr4TWEnTdONQeWO72CXlg32LfctW9i3DzWRjyKVp+hXg48A7WrcbQAR+OUmSv5EkyUiSJN+WJMlfW39MkiQvSJJkguYVqS/TnDfZaVWI3wGemSTJ32k9R2Tz1YWzwCPAS1pD67fS/HBdV6JZXu8+mgsP/+VO/r1JkjyztVjsCcDXgP+1g9jX/RrNTvDvsnmYmyRJXpQkyXjr5pdoDuM/Ycsx061zf3iH5+2GczQ/6P9BkiRXJUkyvj4nOUmSpyVJ8p2tOa5LNP9vi/Tnc+JTNDv+k63329+geZVzay39/wm8qtWhXpkkyeTWJ5KUPfsW+xb7Fj0qTVN/huyH5vD2iQ23n0fzQ/h7Wre/DXgfzQ+FVZof+L+y4fj3A8s0rxTdD3wIGNtw/3mai+02/vynDff/ZOvxCzSv0Hx4Szw/S/MD/yzNjur3gde37isD/xX4JvAVmleAzgPPb93/023OnQIvad3/Bpofit8CvkFz8eLf2OHrt691zj9uc9+/p9kJrNG8ovIe4Motx1xHczHgP9nF/10KvGKb+354m3/7DVuO+7s0P1i/1fr/W/9//0Ga9esfprnw8A+Al+72fdXm/v/eJrb/teH+F7X+bx+kOUz+bppznzc+x1OBz7SO+SZQzfr3yR9//Gn+bP0MwL7FviW1b/EnJWm9yJIkSZLUVU6jkiRJktQTJhuSJEmSesJkQ5IkSVJPmGxIkiRJ6olBrCnsindJ2mwn9ezVnn2LJG3WUd8yiMkG586d29XjSqUSS0tLlz6wj/IWU97iAWPqVN5iyls8MJgxHThwoIvRDLdB6lsgn3HlMSbIZ1zG1Lk8xrXXY9pJ3+I0KkmSJEk9YbIhSZIkqSdMNiRJkiT1hMmGJEmSpJ4w2ZAkSZLUEyYbkiRJknrCZEOSJElST5hsSJIkSeoJkw1JkiRJPWGyIUmSJKknRrIOQHtHvV7g+PECtdp+JicbVCorlMuNrMOSJEkCoFCvM1qtUlhYoDE5yUqlAqVS1mENNZMNdaReLzA9PU6tVgAKAMzPF5mbWzbhkCRJmSvU64xPT1Os1R5tK87Pk3760zA6mmFkw81pVOpItTpKrVbc1FarFalW/eWVJEnZG61WNyUaAMVajcKJE9kEJMBkQx1aWCi0bV9cbN8uSZLUT4WFhbbtyb339jkSbWSyoY5MTrafKjUx4RQqSZKUvcbkZNv29Npr+xyJNjLZUEcqlRWmplY3tU1NrVKprGQUkSRJ0mNWKhVWp6Y2ta1OTdFwGlWmXCCujpTLDebmljl1qkS9vsbEhNWoJElSfjTKZZbn5prVqBYXaUxMsFKpcM3Bg7C0lHV4Q8tkQx0rlxucPt1gaen+rEORJEl6nEa5zPnZ2azD0AZOo5IkSZLUEyYbkiRJknrCaVSSJEnSkGi3y3qjXO7Z+Uw2JEkDLYSwD7gVeCqwD3gd8G3ALwMF4EyM8e3ZRShJ/bHdLuvLc3NQKvXknE6jkiQNtBjjBeCdMcZXAh8E3kIz+bghxvj9wPNDCIezjFGS+mG7XdZHq9WendORDUnSwIsx/kXrrweAZWAVeDCE8EHgauB7gD/Y+JgQwk3ATa3HU9rlVb+RkZFdP7aX8hhXHmOCfMZlTJ3LY1xZxTSyvNy2/arlZehRTJkkG+2GtGOMf97muALwXuA5NIe63xxj/FI/Y5UkDYYQwg8DLwDeSDOJ+AhwM/Ay2oz0xxhvpdlXAaRLu6zTXyqV2O1jeymPceUxJshnXMbUuTzGlVVMY+Pj7GvT/tD4OCNrax3HdODAgY7Pmck0qm2GtNs5ClwRY3w58LPALX0KUZI0QEIIbwD+HvD3YoxfBx6kmXT8MfAqtoxqSNIg2m6X9ZVKpWfnzGwa1ZYh7bPbHHYY+GQI4TXAm4FntTtokIe68xZT3uIBY+pU3mLKWzxgTIMqhPBCmiMUvw98JoTwLZp9yu00p1N9Jsb4hxmGKEl9sd0u6wNbjWrDkPbfvchh7wA+BbwG+LN2BwzyUHfeYspbPGBMncpbTHmLBwYzpp0MdQ+qGOMf0ZyKu9WL+h2LJGWt37usZ1aNasuQ9lqr7WgI4eMbDvsC8MUY47uBFwJf7n+kkiRJknYjk2Rjw5D2d9Ac0v5M665J4NkbDv0o8OQQwueAXwT+374GKkmSJGnXMplGtd2QdoxxFpjdcHsVuKGPoUmSJEnqEjf1kyRJktQTJhuSJEmSesJkQ5IkSVJPZFr6VpIkDbdCvd6s+b+wQGFqisKxYz2t+S+pv0w2JElSJgr1OuPT0xRrtWbDmTOMnznD8tycCYc0IJxGJUmSMjFarT6WaLQUazVGq9WMIpLUbSYbkiQpE4WFhfbti4t9jkRSr5hsSJKkTDQmJ9u3T0z0ORJJvWKyIUmSMrFSqbA6NbWpbXVqipVKJaOIJHWbC8QlSVImGuUyy3NzzWpUi4uMlMssW41KGigmG5IkKTONcpnzs7MAlEolGktLGUckqZtMNiRl4u674eabx1hYKDA52aBSWaFcbmQdliRJ6iKTDUl9V68XuP76ImfPPuHRtvn5InNzyyYckiQNEBeIS+q7anWUs2eTTW21WpFqdTSjiCRJUi+YbEjqu4WFQtv2xcX27ZIkaW8y2ZDUd5OT7adKTUw4hUqSpEFisiGp7yqVFQ4dSje1TU2tUqmsZBSRJEnqBReI01ysWq2Osrw8wvj4mFVxpB4rlxvccccqN9+8xuJigYkJq1FJkjSIhj7ZqNcLTE+PU6sVWy37rIoj9cHBgzA7ez7rMCRJUg8N/TSqanV0Q6LRZFUcSZIk6fINfbJhVRxJkiSpN4Y+2bAqjiRJktQbQ59sVCorTE2tbmqzKo4kSZJ0+YZ+gXi53GBubrlVjeoqxscfsiqOJEmS1AVDn2xAM+GYnT1PqVRiacnqOJIkSXtRoV5ntFqlsLBAY3KSlUoFSqWswxpqJhuSJEna8wr1OuPT0xRrtUfbivPzpJ/+NIxaZTQrQ79mQ5IkSXvfaLW6KdEAKNZqFE6cyCYgASYbkiRJGgCFhYW27cm99/Y5Em3kNCpJkqSc27oWgZMnnRq0RWNysm17eu21fY5EG2WSbIQQrgJ+AjgOnIgxfvgixz4E3Nm6eSrGeHvvI5QkScqHdmsR0rvuonDbbTTK5Qwjy5eVSoXi/Pym12l1aorUaVSZympkYwL4JnBbB8cuxBhfcbEDQgg3ATcBxBgp7bLqwMjIyK4f2yt5iylv8YAxdSpvMeUtHjAmSfnUbi1CcvYso9Uq52dnM4oqfxrlMstzc80RoMVFGhMTrFQqXHPwICwtZR3e0Mok2Ygx1oAPhBBOdHD4Wgjh94Al4O0xxrNtnu9W4NbWzXRpl2+oZunbfL0Z8xZT3uIBY+pU3mLKWzwwmDEdOHCgi9FIysJ2axEKi4t9jiT/GuWyCVjO5H6BeIzxmTHGlwEfAX4p63gkSZL6abu1CI2JiT5HIu1crpKNEMLREMLHt7k7Af6yn/Fos3q9wI03Fjh6dD8zM2PU64WsQ5IkaeCtVCqsTk1taksPHWpuWCflXFYLxJ8GfAI4ADwcQviBGOPrgEng2RuOewrwm8DDwH3AWzMIVzQTjenpcWq1AtBMMubni8zNLVMuN7INTpKkAdZuLcLIyZM0rEalPSCrNRv3AC9s0z4LzG64fR/w0j6Gpm1Uq6PUasVNbbVakWp1lNnZ8xlFJUnSkEnTrCOQdsR9NtSRhYX2U6YWF51KJUlSL1n6VntZrtZsKL8mJ9tPlZqYcAqVJEm9dLHSt1LemWyoI5XKClNTq5vapqZWqVRWMopIkqThYOlb7WVOo1JHyuUGc3PLnDpVol5fY2KiQaWy4uJwqYfq9QLHjxeo1fYzOenvnDSsLH2rvcxkQx0rlxucPt1gaen+rEPJNb8gqhusACdp3UqlQnF+fvOaDUvfao8w2ZC6yC+I6hYrwElaZ+lb7WUmG1IX+QVR3WIFOA2LQr3e/BK9sEBhaorCsWNWWGqjUS5zfvbR3QEolUqwtJRhRFJnTDbUMacHXZpfENUtVoDTMHhcSdczZxg/c4bluTkTDmlAWI1KHVmfHjQ3V+DMmSu5/fZ9TE+PU6/7JXojvyCqW6wAp2HQrqRrsVazpKs0QEw21JGLTQ/SY/yCqG5ZrwA3Pd3g8OGHOXLkgmt/NHAs6SoNPqdRqSNOD+qMJYLVTVaA06CzpKs0+Ew2aE4RqlZHWV4eYXx8zC+HbTg9qHN+QZSkzrQr6bo6NWVJV2mADH2y8Vip0vUpQvssVdpGpbLC/Hxx01Qqpwe150J6SepMo1zm/Hvew9ixYxQeeACuuYbz73mPi8OlATL0yYalSjvj9KDOuM+GJHWuUK8z9ra3UbznnmbDAw8w9ra3WY1KGiBDv0DctQg7l6ZZR5BfLqSXpM5ZjUoafEM/suFahM54xb4zJq+S1DmrUUmDb+hHNiqVFQ4c2Fyq9MAB1yJs5RX7zpi8SlpXqNcZm5lh5LrrGJuZoVCvZx1S7liNShp8Qz+yAZAkyUVvyyv2nXIhvZQ/IYSrgJ8AjgMnYowfDiF8L/DLNIdqz8QY397Nc27dGXsfUJyfdy3CFlaj6lyhXme0WqWwsNBM0k6ehFEv+Cn/hn5ko1od5Wtf25xzfe1rI16x38Ir9p1xI7bO3X03zMyMcfTofmZmxtyNXr00AXwTuG1D263ADTHG7weeH0I43M0TuhahM+vVqFaf9jQeefKTeWRqympUbawnr/tuv50rz5xh3+23U3zVqxwt054w9CMbXrHvjFfsO+c+G5dWrxe4/voiZ88+4dE21wCpV2KMNeADIYQTACGEa4BV4MEQwgeBq4HvAf5g4+NCCDcBN7Weg1Kp1PE5R5aX27Zftby8o+fppZGRkexjuftuij/zMyQbqlGVfuZnWL3jDjh4MNvYNsj6tSocP05hS/KanD1L6dQpGqdPZxTV42X9Om0nj3ENU0xDn2x4xb4zlr5VN1Wro5w9u3m6oiWn1WcHgI8ANwMvo81If4zxVpojIADp0tJSx08+Nj7OvjbtD42Pc34Hz9NLpVKJnfybemHs5pt5wtmzm9qSs2dZu/lmzs/OZhTV42X9Wu2v1Wh3CXStXuf+nLyfIPvXCR4/3WylUuGav/W3Mo9rqzy8VlvtJKYDBw50/LxDP42qUllhamrzAnGv2F+cpW91uRxRVJZijH8JPAi8Efhj4FVsGdW4XCuVCqtTU5vaXIvweFaj6owL6TvTbrrZ+PR0c96uMjP0IxvrV+yr1VGWl69ifPwhr9i3YelbdZMjiuqnEMLTgE/QHM14OITwA8CbgdtpTqf6TIzxD7t5zka5zPLcHKPVKlctL/PQ+DgrlYprEbbwS3RnvnHDDTzxt3+bZG3t0bZ0ZIRv3HBDhlHlz3ZrpRonTsAtt2QTlEw2oJlwzM6ebw0fOYWjHXdaVzdVKivcddcTN02lckRRvRJjvAd4YZu7XtTL8zbKZc7PzlIqlXIzdSpv/BLdmSd95CObXiOAZG2NJ33kI5x/8Ysziip/thspS+69t8+RaKOhn0alzjjtRd1ULje4445Vjhy5YNUuaYhd7Eu0HuN0s85sN1KWXnttnyPRRo5sqCNOe1G3HTyIo2LSkPNLdGecbtaZ7fZtSU+cyC4oObKhzriQXpLUbX6J7ky7ggPpoUMWHNhifa3UhSNHePjwYS4cOcLy3FyuyigPI0c21BFL30qSus0dxDuzseBAYXGRxsQEIydP0nAH8cdZXyul/Mgk2QghXAX8BHAcOBFj/PA2xxWA9wLPoVkC6c0xxi/1K05t5mZ1kqRu2voleqRcZvnYMat2tbH1S3SpVAILD2gPyGpkYwL4JnDbJY47ClwRY3x5COElwC3Adb0OTpIk9cfGL9GlUomGX6ClgZJJshFjrAEfCCGcuMShh4FPhhBeQ7Mm+rPaHRRCuAm4qfXcu95qfZi2jt+tvMUDxtSpvMWUt3jAmCRJ6ra9sGbjHcCngNcAf9bugBjjrcCtrZvpbrd/3+tbx/dD3uIBY+pU3mLKWzwwmDEdOHCgi9FIkrQzuapGFUI4GkL4+IamLwBfjDG+m+aGTF/OJjJJknamUK8zNjPDyHXXMTYzQ6FezzokaeCt/97tP3rU37ucyGqB+NOATwAHgIdDCD8QY3wdMAk8e8OhHwWuCyF8DlgDfrLvwUo7VK8XOH68QK22n8lJq3ZJw6hQrzM+Pf1olaV9QHF+nuW5ORc/Sz2y9fcOmr936ac/DVbuykxWazbuoTlSsbV9FpjdcHsVuKGPoUmXpV4vMD09Tq1WoFlADebni+6OLQ2Z0Wp10xcegGKtxmi1allO7UqhXm9W7VpYaO5PcvKkX6C32O73rnHiBNxySzZBaU+s2ZD2jGp1lFqtuKmtVitSrY66W7Y0RNwZW93U7op9etddFG67zZGyDbb7vUvuvbfPkWijXK3ZkPa6hYVC2/bFxfbtkgbTI9tccX7k6qv7HEn+bZxjX7jxRufYt9Huin1y9iyj1WpGEeXTdjvSp9de2+dItJEjG1IXjY4+0rb96qvbt0vSMHvcFfszZxg/c8a1LVs4UtaZ7XakT0+cyC4oObIhSVK3XbGy0r79wQf7HEm+XWxtix6z3RX7xsREnyPJt/Ud6S8cOcLDhw9z4cgRlufm4ODBrEMbao5sSF20stI+f3/wQfN6aZj45bAzXrHvTLsr9umhQ6xUKhlGlU8bd6RXPvgNSB2r1wvceGOBo0f3MzMzRr3uOoStJifbV5yamLASlTRMVioVVqemNrWtTk355XAL17Z0pt0V+9U77nCqmfYERzbUEUu6dqZSWWF+vripItXU1CqVSvspFZIG0/qXw9FqlauWl3lofJyVSsUvh9q1rVfsS6USLC1lGJHUGZMNdcSSrp0plxvMzS1z6lSJen2NiQk39ZOG1fqXw1KpxHm/FLbl2hZp8JlsqCOWdO1cudzg9OkGS0v3Zx2KJOWaa1ukweeaDZpThGZmxrjuuhHXImzDtQiSpG5zbYs0+IZ+ZOOxtQjrU4T2uRahjUplhc9/vsi5c49NpTpwwLUIkqTd27i2pbC4yEi5zPKxY65tkQbI0CcbrkXoXJIkF70tSdJObVz4XCqVaLi+RRooQz+NyrUInalWR/na1zbnpl/72gjVavuyhZIkSdLQj2y4FqEzJmXqtrvvhptvHmNhocDkpFW7JEkaREOfbLgvQmdMytRN9XqB668vcvbsEx5tc62UBk2hXme0WmVkeZkx99mQNKSGPtlY3xehWh1lefkqxscf8gprGyZl6qZqdZSzZzev+XGtlAZJoV5nfHqaYq0GwD6gOD/P8tycCYd2ZT15LSwsNEsGnzwJ2+zALuXJ0Ccb0Ew4ZmfPUyqVWFryi047blbXuXq9wPHjBWq1/U4P2obT8jToRqvVRxONdcVajdFqddMu0FIntiavAOldd1G47TaTV+WeyYY65mZ1l/ZYKeUC0Pzi7PSgx3NangZdYWGhffviYp8j0SBol7wmZ8+avGpPGPpqVFI3XayUsh5Tqaxw6FC6qc1peRok7oytbjJ51V5msiF1kdODOlMuN7jjjlWOHLnA4cMPc+TIBUd/NFDcGVvdZPKqvcxpVFIXOT2ocwcP4mJwDayNO2NftbzMQ1aj0mVYqVQozs9vXrNx6JDJq/YEkw2pi6zaJWnd+s7YpVKJ8+6KrcuwMXktLC7SmJhg5ORJGlaj0h5gsiF1kVW7JEm9sJ68riuVSmASqz3AZEPqMqt2SZK6zX02tFeZbEiSJOWY+2xoL7MalaRM3H03zMyMcfTofmZmxqjXrdglSe1cbJ8NKe8c2ZDUd/V6geuvL3L27BMebXPzQ0lqz302tJc5siGp76rVUc6eTTa1ufmhJLXnPhvay0w2JPWdmx9qGBTqdcZmZhi57jrGZmYo1OtZh5RL66/T/qNHKdx4o69TG+02iXSfDe3Wxt+5fnw2ZTKNKoRQAN4LPAcoAG+OMX5pm2MfAu5s3TwVY7y9P1FK6hU3P9Sg27qgdx9QnJ9neW7OBb0bPG7h85kzjJ854+u0hftsqFvaFRtY/2yiVOrJOXc9shFCuPPSR23rKHBFjPHlwM8Ct1zk2IUY4ytaPyYa0gCoVFY4dCjd1Obmh7qYy+xz+q7dgt5ireaC3i18nTq3vs/G/R//eHO/jYMHsw4pl/p91X6vyeJ37qIjGyGEv36Ru9tPIOzMYeCTIYTXAG8GnnWRY9dCCL8HLAFvjzGebRPnTcBNADHG5kY3uzAyMrLrx/ZK3mLKWzxgTJ3KU0ylEnzmMynveMcj3HtvwrXXppw4kXLw4DWZxpWn12hdHmPqlR72OX3ngt7O+Dqpm7a7ap9++tPuSdKSxe/cpaZR/XfgHiBpc9+1l3nudwCfAl4D/Nl2B8UYnwkQQvgR4JeAv9/mmFuBW1s306Vd7qhZKpXY7WN7JW8x5S0eMKZO5S2m7/iOErfcsjmerMPL22sElx/TgQMHuhhNz/Wyz+krF/R2xtdJ3bTdVfvGiRNwy8Um0QyPi/3O9WptxaWetxZj/Gvt7gghfOUyzvsFoBhjfHcI4fuAL7ee8yjw2hjj32vzmAT4y8s4pyQp33rV5/TdSqVCcX5+0xef1akpF/Ru4eukbtruqn1y7719jiS/LvY716u5BZdKNuJF7vsfl3HejwLXhRA+B6wBP9lqnwSevX5QCOEpwG8CDwP3AW+9jHNKkvKtV31O321c0HvV8jIPjY+zUqm46HmLrQufR8pllo8d83XSrmx31T69dk8NjPZUu2IDvf5sStI0vfRRe0t67ty5HT2gXi9QrY6yvHwV4+MPUams5GZjsbxN68hbPJC/mOr1AqdOlajV1picbOTm/ZS31ynD42KVAAAgAElEQVRv8cBgxtSaRtVuWpJ2Zsd9y7o8vq8gn3HlMSbIZ1zG9HiFep39R48y8rWvPdq29tSn8sjv/A5LOVuzkfVr1c5OYtpJ39Lx9KwQwpOA/RvbYox7fol/vV5genqcWq3YatnnTsbatcfeTwWaVZ3dGVvajUHtcyT11taL6AN4UX3P6aj0bQjhp4BzwP8CPg98Bbijh3H1TbU6uiHRaHInY+2W7yfp8g1ynyOpd0arVYpbRiCL585ROHEim4AEdL7Pxk8BU8BfxBivBd4I/H7PouojdzJWN/l+krpiYPscSb3jAvF86jTZSGKM54FHAGKMvw68smdR9ZE7GaubfD9JXTGwfY6k3nGBeD51umZjJYTwBOCuEMKvAv+zhzH1VaWywvx8cdPUF3cy1m75fpK6YmD7HEm9s11Z19RpVJnqdGTjeporzn+aZoLySuD1PYqpr8rlBnNzyxw5coGXv/wRjhy54GJe7dr6+2l6usHhww/7fpJ2Z2D7HEm9s17W9cKRIzx8+DAXjhxheW4ODh7MOrRcKdTrjM3MsP/oUcZmZijUe1t7o6ORjRjjl1t/vQ94U+/CyUa53GB29nyr5Nf5rMPRHlcuNzh9usHS0v1ZhyLtSYPe50jqnUa5zPnZ2azDyK1Cvc749PSm0Z/i/HwzKSuVenLOjpKNEMKHgMfVDosx/njXI5IkDbVB6XMK9Tqj1Sojy8uMuamfpBwYrVY3JRoAxVqN0WoV5uZ6cs5O12z80Zbbrwc+191QlHf1eoHjxwvUavtztVld3vg6SZdtz/c5W68e7uOxq4cmHJKyUvjqV9u3b0lAuqnTaVS/uvF2COHXgdt7EpFyyc3qOuPrJF2+fvQ5IYQC8F7gOTR/Wd8cY/xSt57/YlcPneIhKStXfP3r7dvvu+/xw8ndOucuH5cC7euLaSC5WV1nfJ2knuhFn3MUuCLG+HLgZ4Fbuvnk29X7LywudvM0krQj6VOesqP2buh0zcbX2Tx/dhT45Z5EpFxys7rO+DpJl69Pfc5h4JMhhNcAbwae1SaOm4CbAGKMlHaweLIwNQVnzjyufaRc3tHz9NLIyEhuYlmXx5ggn3EZU+fyGFdWMRW+8zthfr5te9KjmDpds/HCDX9PgftjjN/oejTKLTer64yvk9QV/epz3gF8CngN8Gdb74wx3grcuh7H0tJSx09cOHaM8TNnHlfvf/nYMRo7eJ5ealZgzEcs6/IYE+QzLmPqXB7jyiqmi302XbO21nFMBw4c6Picna7Z6N2qEe0JblbXGV8n6fL1qc/5AlCMMb47hPB9wJcv9YCdWK/3P1qtctXyMg9ZjUpSDmz8bCosLtKYmOj5Z9NFk40QwiO0KT/YaluJMV7Tk6iUO+ub1Z06VaJeX2NiwipL7fg6SbvX5z7no8B1IYTPAWvAT3bxuYHH6v2XSiXO5+yqqqTh1e+9SC41sjFKcxfXPwWeCzSAIs0KHv+xt6Epb9ysrjO+TtKu9a3PiTGuAjd08zklSY930WRjfY5sCOFCjPGB9fYQwhuAu4Bf7214kqRhYZ8jSYOn0wXi/y2E8C7g/cDDwHU0rzZJktRt9jmS1COFer25ZmNhgcbkZLZrNjY4RrMG+TzNjVDvwuFnSVJv2OdIUg8U6nXGp6c3VaMqzs+zPDcHPSrF22k1qvPAG3oSgSRJG9jnSFJvjFarmxINgGKtxmi1CnNzPTnnbncQJ4RQ6WYgkobL3XfDzMwYR4/uZ2ZmjHrdjQ+1PfscSbp8hYWF9u2Liz07Z6fTqNp5I1DtViCShke9XuD664ucPfuER9vm54vMzS1bJniDer3A8eMFarX9TE4OfRll+xxJukyNycn27RMTl5UUXMyl9tn4DeCrwMu23JUAT+tRTJIGXLU6ytmzyaa2Wq1ItTrK7Oz5jKLKl3q9wPT0OLVaAWiO+gx6QmafI0m9tVKpUJyff9wO4iuVCr3aPO9SScw54D7gWjYvzkuwBKGkXVpYaD9lanHRqVTrqtXRTTvRw1AkZPY50ja2VhDi5EkYHc06LO0xudtBPMb4/wCEEH4sxvi5jfeFEB7qWVSSBtrkZPsr8xMTg3nFfjeGMSGzz5Haa1dBKL3rLgq33dbTL4kaTP3eQbyjBeIxxme3af6hLsciaUhUKit8x3ekm9qe+tQ1KpWVjCLKn2FOyAalzyneeSff/qIXUXzKU5p/3nln1iFpj2pXQSg5e7ZZQUjKuV1Xo4ox/lk3A5E0XNJ06+20/YFDqlJZYWpqdVPb1NTq0CZke63PKd55J6XXvpbiPfeQ/NVfUbznnuZtEw7tQhYVhKRu6SjZCCG8OoRwRevvbw4h/NsQwnN6G1r/1OsFZmbGuO66EUtwSn1QrY5yzz2bF4ifO9dcj6CmcrnB3Nwy09MNDh9+mCNHLgz04vCNBqHPGTt2jGRtbVNbsrbG2LFjGUWkvexiFYSkvOu0ytWvxBg/GUL4TqAC/CLwPuClPYusTx6r+LK+EHPfwFd8kbI2jOsRdqNcbnD6dIOlpfuzDqXf9nyfU3jggR21SxfTroJQeugQKxW3n1H+dZpsfLP1548C/zLG+MEQws27PWkIoQC8F3gOzZqOb44xfmm3x12OIa34ImVqmNcjqCNd7XOy0Hjyk7miTWLRePKTM4hGe127CkIjJ0/SsBqV9oBOk43FEMLbgB/nsStLu17vARwFrogxvjyE8BLgFuC63R4XQrgJuAkgxkipVOo4kOXl9i/B8vJVO3qeXhkZGclFHOvyFg8YU6fyFNPJk3DXXemmvTYOHUo5eTLbGPP0Gq3LY0x90O0+p+/OnzpF6bWv3TSVKh0Z4fypUxlGpb1sawWhUqkES0sZRiR1ptNk403AzcDbY4y1EMLVwH+5jPMeBj4ZQngN8GbgWZdzXIzxVuDW1s10aQe/fOPjY8C+Nu0PsbSU/chGqVRiJ/+eXstbPGBMncpTTKOjcMcdJW6+eY3FxQITE83dsUdHG5n2nXl6jdZdbkwHDhzoYjR90+0+p+9WX/xilj72McaOHWNkZYW10VHOnzrF6otfnHVoktRXHSUbMcb/CfzYhtsPAjde5rnfAXwKeA1wsSojnR63K5XKCvPzxU1TqYa54ovULwcP4lTFS6jXCxw/XqBW28/kZDMhG4a1ZD3qc/pu9cUv5uuf/3wuk1hJ6peOh6VDCC8NIfxY6++FEMK3X8Z5vwB8Mcb4buCFwJdbz3s0hPDxSx3XTesVX44cucDLX/7IUFV8kZRf68Ur5uYKnDlzJbffvo/p6fGhqZbX5T5HkpSRjkY2Qgg/D/wdYBL4EJAAvwM8b5fn/ShwXQjhc8Aa8JOt9kng2R0c11XlcoPZ2fOtq09eaZWUvWEuXtGDPkeSlJFO12y8DvibwBcBYoxrIYTixR+yvRjjKnBDm/ZZYPZSx0nSoPvqV9uPYNRqQzGy0dU+R5KUnU6nUX0rxvjQ+o0QwihwVW9CkiR9/evtP57vu29PFWXaLfscSRoQnfZanwwhnAKeFEL4B8B/AH67d2FJ0nB7ylPSHbUPGPscSRoQnSYbPw/cBfw34EeA3wTe3qugJGnYTU2t7ah9wNjnSNKA6LT07SPAv2r9ABBCSLZ/hCTpcgxzWW77HEkaHBdNNkIIBeAfAE8H/kuM8Xda7dPAO9lcOUqS1CXrZblPnSpRr689uvHhIJflts+RpMFzqZGNW4FnAncC/zyE8GvANPAE4G09jk2SBKRDsUwDsM+RdJkK9Tqj1SqFhQUak5OsVCpQKmUd1lC7VLLxMuCvxxhXQwj/BFgE3hJj/EDvQ5Ok4bW+qV+z1G2z3O38fHHQNx0dqD5n/UvPyPIyY+PjrFQqNMrlrMOSBlahXmd8eppirfZoW3F+nvTTn4bR0QwjG26dLBAfCSHsAx4GasBtIYR9rTZJUg9cbFO/ATcQfc76l559t9/OFZ/7HPtuv53x6WkK9XrWoUkDa7Ra3ZRoABRrNQonTmQTkIBLj2w8A3iQ5u6t69Zvp6xfbpMkddXCQvuP18XFgf7YHZg+Z7svPaPVKudnZ7d5lLS9rdODOHnSq/VbFBYW2rYn997b50i00UWTjRjjUOweJUl5MznZfqrUxMTATqEaqD5nZEuical26WLaTQ9K77qLwm23OTVvg8bkZNv29Npr+xyJNhqYD3ZJGiSVygpTU6ub2oal9O0gSO67b0ft0sW0GylLzp5ltFrNKKJ8WqlUWJ2a2tS2OjVFw2lUmeponw1JUn8NY+nbQfLIt3873HPP49uf8pQMotFet930oMLiYp8jybdGuczy3FxzutniIo2JCVYqFa45eBCWlrIOb2iZbEhdVq8XOH68QK22n8lJvyBq98rlBqdPN1hauj/rULRDjac/Hb74xce3b7nqKnViu+lBjYmJPkeyhwxRzfC8M9mQumhIy5VK2mKlUqE4P79p6svq1FSz5r+0Q+3eT+mhQ76ftrD0bT65ZkPqoiEuVyppg/XpHBeOHOGRl7+cC0eOsDw352Je7crG99PDhw9z4cgRVu+4w/fTFpa+zSdHNqQuGtJypZLaaJTLnJ+dpVQqcd754rpM6++ndaVSyXUIW1j6Np8c2ZC6aBjLlUqSlAePbDNVKnUKVaZMNqQuslypJEnSY5xGJXWR5UolScrGFSvtL+wl27SrP0w2pC6zXKkkSf3nDuL55DQqSZIk7XnuIJ5PjmxIkiRpz3MH8Xwy2ZAkSdJA2FoiWNlzGpUkSZKknnBkQ5IkSQOhUK83p1EtLNCYnGSlUoFSKeuwhprJBnDnnUWOHRtjZWWE0dFv59Sp87z4xauXfqAkSVIfbP0SzcmT4GZ1mxTqdcanpynWao+2FefnST/9aV+rDA39NKo77yzy2teWuOeeIn/1Vwn33NO8feedxaxDkyRJevRL9L7bb+fKM2fYd/vtFF/1Kgr1etah5cpotbop0QAo1moUrEaVqaFPNo4dG2NtLdnUtraWcOzYWEYRSZIkPabdl+jk7FlGq9WMIsqnwsJC2/bk3nv7HIk26vs0qhDCq4GbgQLwsRjjP7vIsR8Gng+cB/48xvimbsfzwAOFHbVL6o6774abbx5jYaHA5KQ7rUvSdrb7El1YXOxzJPnmpn751NdkI4SwD6gCLwa+AfxxCOE3Y4wXGwf86RjjZ3sV05Of3OCBBx4/wPPkJ/ulR+qVer3A9dcXOXv2CY+2zc8XmZtbNuGQpC22+xLdmJjocyT5tlKpUJyf3zQKtDo1Reo0qkz1LNkIIbwf+O4tzceBPwW+Dfgg8DDwAmC7ZGMJ+KchhG8C748xfnSbc90E3AQQY6S0g6oDH/pQyg/+YLppKtXISMqHPpTu6Hl6ZWRkJBdxrMtbPGBMncpTTMePFzh7dvP0xVqtyKlTJU6fzi7ZyNNrtC6PMUnqr3ZfotNDh5qVlvQoN/XLpyRN076dLITwAiACX6KZeMwAvxdj/K1LPG4MmAeeG2P8xiVOk547d25HcW2uRrWWq2pUpVKJpRz9guQtHjCmTuUppqNH93PmzJWPaz98+GE+/vH7M4ioKU+v0brLjenAgQMAyaWO0yXtuG9Zl8f3FeQzrjzGBPmI69FqVK0v0SMnT7KUswpLeXid2sljXHs9pp30Lf1es/EVmovSbwDWgFcAtwCEEK4BPgG8Ncb4J1se16A5CvKtXgT14hev8vnPfz2X//HSIJqcbD96MTHhFCpJamfrztilUsmr9doT+ppsxBgfDCG8A/gssAq8L8a4fqnoicB3AY+WgQohnAKeSzPZeGuMMR/DDZIuS6Wywl13PXHTVKqpqVUqlZUMo5IkSd3W92pUrXUXj1t70Uo6Jra0HetXXJL6p1xu8P73r/JjP5bwwAMFnvzkBu95z3kXh0uSNGCGfp8NSf1Xrxd405uK3HNPkQceuIJ77inytreNUa9bclqSpF4q1OuMzcyw/+hRxmZmer45ZN9HNqRBV68XOH68QK223/0jtlGtjratRlWtjjI7ez6jqDSIQggvA34OuDLG+MpW25OAXwMO0Jyme2OM8Z7sopSk/ljfjX5jZbPi/DzLc3PQo8qHjmxIXVSvF5ieHmdursCZM1dy++37mJ4e94r9FgsL7V+PxUVfJ3Xd84BfYHPVlBngrhjjK4APAf8og7gkqe/a7UZfrNV6uhu9IxtSF1Wro9RqxU1tXrF/vNHRR9q2X311+3bpUrbZ2+kNMcbZEMLTt7QfBn4+hPAG4IeBp23znLvew2mjvO6Vkse48hgT5DMuY+pcHuPKKqaR5eW27VctL0OPYjLZkLrIK/ZSNmKMb9rhQ2aBjwHXA5/b5jlvBW5t3Ux3Wxo9r2XV8xhXHmOCfMZlTJ3LY1xZxTQ2Ps6+Nu0PjY8zsra20302OuI0KqmL3D+iMysr7T96HnzQjyT1xReAO2KM7wP+b+D3M45HkvpipVJhdWpqU9vq1FRPd6N3ZEPqokplhfn54qapVO4f8XgmZeqXEMIHgO8FDoYQ/gh4LXAK+EgI4VXAg8CPZxiiJPVNo1xmeW5u0270K5UKjXK5Z+c02ZC6qFxuMDe3zKlTJer1NSYmrEbVjpv6qV9ijG/c5q4f6msgkpQTW3ej7zWTDanLyuUGp083WFq6P+tQcqtcbnDHHavcfPMai4sFkzJJkgaUyYakTBw8iBW6JEkacK7GlCRJktQTjmxIysTdd8PNN4+xsFBwp3VJkgaUyYakvqvXC1x/fZGzZ5/waNv8fJG5uWUTDkmSBojTqCT1XbU6uqkSFTy207qk4VKo1xmbmWH/0aMUbryRQr2edUiSusiRDUl9507rnanXCxw/XqBW2+9UMw2kQr3O+PQ0xVqt2XDmDONnzrA8N9fTuv+S+sdkQ1LfuanfpdXrBaanx6nVCkAzCXOqmQbNaLX6WKLRUqzVGK1W+7oPgKTecRqVpL6rVFY4dCjd1OamfptVq6ObdqIHp5pp8BQWFtq3Ly72ORJJveLIhqS+c1O/S3OqmYZBY3KyffvERJ8jkdQrJhuSMuGmfhfnVDMNg5VKheL8/KapVKtTU6xUKhlGJambnEYlSTlUqawwNbW6qc2pZho0jXKZ5bk5Lhw5wsOHD9OYnnZxuDRgHNmQpBwqlxvMzS1z6lSJen3NqWYaWI1y+dHF4KVSicbSUsYRSeomkw2aVV+q1VGWl0cYHx+zQ9dlsVypuqVcbnD6dIOlpfuzDkWSpF0Z+mTjsfKS61Vf9lleUrtmuVJJkqTHDP2aDctLqpt8P0mSJD1m6JMNy0uqm3w/SZIkPWbop1FZXlLd5PtJknamUK8zWq1SWFigMDVF4dgxq1FJA2ToRzYqlRWe+tS1TW1Pfeqa5SW1K5YrlaTOFep1xqen2Xf77Vx55gyFuTnGp6cp1OtZhyapS4Y+2QBI0/Sit6VOrZcrnZ5ucPjwwxw5csHF4ZK0jdFqddOGfgDFWo3RajWjiCR129BPo6pWRzl3bvOC3nPnmgt63d1Yu2G5UknqTGFhoX374mKfI5HUK31PNkIILwN+DrgyxvjKSxz7auBmmjVEPxZj/GfdjscFvZIkZaMxOdm+fWKiz5FI6pUsplE9D/gFILnYQSGEfUAV+EHg+4E3hhC6vmLMBb2SJGVjpVJhdWpqU9vq1BQrlUpGEUnqtp6NbIQQ3g9895bmN8QYZ0MIT+/gKZ4F/CnwbcAHgYeBFwCPWzUWQrgJuAkgxkipVOo4zpMn4a67Us6efSz3OXQo5eTJkR09T6+MjOQjjnV5iweMqVN5iylv8YAxSf3WKJdZnptrVqNaXGSkXGbZalRST22sANeYnGSlUunp71zPko0Y45u68DTPB04Bx4GZi5zrVuDW1s10aWmp4xOMjsJttxWoVkdZXr6K8fGHqFRWGB1tsIOn6ZlSqcRO/j29lrd4wJg6lbeY8hYPDGZMBw4c6GI0Uvc1ymXOz84Czfd7I2e/g9IgWa8At7EwQ3F+nuW5OejRha3cLBAPIVwDfAJ4a4zxT4Cv0JzmdQOwBrwCuKUX5y6XG8zOnm916i4KlyRJ0uC5aAW4ubmenLPvazZCCB8A/h3w3SGEPwohPKN11xOB7wLGAGKMDwLvAD4L/GfgfTHGc/2OV5IkSRoEWVSA6/vIRozxjdu0nwMmtrR9FPhoP+KSJEmSBtnFKsD1KilwUz9JkiRpCGRRAS43azYkSZIk9c7WCnCNiYm9W41KkiRJUr5srADXD06jkqScqtcL3HhjgaNH9zMzM0a9Xsg6JEmSdsSRDUnKoXq9wPT0OLVaAWgmGfPzRebmlimXG9kGJ0lShxzZkKQcqlZHqdWKm9pqtSLV6mhGEUmSBkGhXmdsZob9R48yNjNDoV7v6fkc2ZCkHFpYaD9lanHRqVSSpN3JYgdxRzYkKYcmJ9tPlZqYcAqVJGl3LrqDeI+YbEhSDlUqK0xNrW5qm5papVJZySgiSdJeNxQ7iEuSLq1cbjA3t8ypUyXq9TUmJhpUKisuDpck7VoWO4ibbEhSTpXLDU6fbrC0dH/WoUiSBsBKpUJxfn7TVKr1HcSv6dE5TTYkSZKkIeAO4pKkR9XrBY4fL1Cr7Wdy0mlUkqTL1+8dxE02JCmH3NRPkjQIrEYlSTnkpn6SpEFgsiFJOeSmfpKkQWCyIUk55KZ+kqRBYLIhSTnkpn6SpEHgAnFJyiE39ZMkDQKTDUnKKTf1kyTtdU6jkiRJktQTJhuSJEmSesJkQ5IkSVJPmGxIkiRJ6gmTDUmSJEk9YbIhSZIkqSdMNiRJkiT1hMmGJEmSpJ5wUz9J0sAKIbwRuB4YB94fY/wXIYQnAb8GHAAawI0xxnsyDFOSBlbfk40QwsuAnwOujDG+8hLHfhh4PnAe+PMY45t6H6EkaYD8LvBB4GrgL4B/AcwAd8UYfzSEcAPwj4CbsgtRkgZXFiMbzwN+AfilDo//6RjjZy92QAjhJlodRYyRUqm0q8BGRkZ2/dheyVtMeYsHjKlTeYspb/GAMe1lIYT3A9+9pfkNMca7Wvc/Bai32g8DPx9CeAPww8DTtnnOge1bIJ9x5TEmyGdcxtS5PMY1TDH1LNm4yAf/bAjh6R0+zRLwT0MI36Q5/P3RdgfFGG8Fbm3dTJeWlnYTMqVSid0+tlfyFFO9XuDUqRK1WsrkZINKZYVyuZF1WLl6jdYZ06XlLR4YzJgOHDjQxWjy62Ij3yGECeBfAT++oXkW+BjNKVaf2+Y5B7ZvgXzGlceYIJ9xGVPn8hjXXo9pJ31Lz5KNbkx5ijG+HSCEMAbMhxD+XYzxG5cdnHasXi8wPT1OrVYACgDMzxeZm1vORcIhSe2EEA7RTDTeEmP801bzF4CHYozvCyH8feD3MwtQkgZcbhaIhxCuAT4BvDXG+Cdb7m4ADwPf6ntgAqBaHaVWK25qq9WKVKujzM6ezygqSbqk36LZ1/1qCAHg3cAp4CMhhFcBD7J5xEOS1EVZLBD/APC9wMEQwh8Br40x/gXwROC7gLENx54Cnksz2XhrjHG13/GqaWGh0LZ9cbF9uyTlQYzxudvc9UN9DUSScqJQrzNarVJYWKAxOclKpUKjXO7Z+fqebMQY37hN+zlgYkvbsb4EpUuanGw/VWpiwilUUq/U6wWOHy9Qq+3P1TopSdLeVKjXGZ+eplirPdpWnJ9neW4OerRgPTfTqJRvlcoK8/PFTVOppqZWqVRWMoxKGlyuk5Ikddtotbop0QAo1mqMVqswN9eTc7qDuDpSLjeYm1tmerrB4cMPc+TIBb/0SD10sXVSkiTtRmFhoX374mLPzunIhjpWLjc4fbrB0tL9WYciDTzXSUmSuq0xOdm+fWKiZ0mBIxuSlEOuk5IkddtKpcLq1NSmttWpKVYqlZ6d02RDknKoUllhampzAT7XSUmSLkejXGZ5bo4LR47w8OHDXDhyhOW5ucGqRiVJurT1dVKnTpWo19eYmLAalSTp8jXKZc7PzvbtfCYbkpRTrpOSJO11TqOSJEmS1BMmG5IkSZJ6wmRDkiRJUk+YbEiSJEnqCZMNScqper3AjTcWOHp0PzMzY9TrbugnSdpbrEaljtXrBY4fL1Cr7Wdy0jKcUi/V6wWmp8ep1QpAM8mYny8yN7fs750kac8w2VBH/OIj9Ve1OkqtVtzUVqsVqVZHmZ09n1FUkqS9rlCvM1qtUlhYoDE5yUql4qZ+yp5ffKT+WlhoP2VqcdGpVJKk3SnU64xPT1Os1R5tK87Pszw3B6VST87pmg11xC8+Un9NTrYfMZyYcCRRkrQ7o9XqpkQDoFirMVqt9uycJhvqiF98pP6qVFaYmlrd1DY1tUqlspJRRJKkva6wsNC+fXGxZ+c02VBH/OIj9Ve53GBubpnp6QaHDz/MkSMXXCMlSbosjcnJ9u0TEz07p2s21JH1Lz6nTpWo19eYmLAaldRr5XKD06cbLC3dn3UokqQBsFKpUJyf3zSVanVqipVKhWt6dE6TDXXMLz6SJEl7V6NcZnlurlmNanGRxsSE1agkSZIkdUejXOb87GzfzueaDUmSJEk9YbIhSZIkqSdMNiRJkiT1hMmGJEmSpJ4w2ZAkSZLUEyYbkiRJknrCZEOSJElST/R9n40QwhuB64Fx4P0xxn9xkWNfDdwMFICPxRj/WX+ilCRJknS5shjZ+F3gB4CXACe2OyiEsA+oAj8IfD/wxhBC77Y3lCRJktRVSZqmPXniEML7ge/e0vyGGONdrfufQXO04oXbPP4FwM8BbwPeAzwDeFeM8bfaHHsTcBNAjHHrOSVp2CVZBzAAetNZStLe1VHf0rORjRjjm2KML9zys55oTAD/CvjxSzzN84FTQAX4zxc5163r56D5D9/VTwjhC5fz+F785C2mvMVjTHs3przFM+Ax6fJl/X+Y1/fWwMeU17iMaW/HNSAxdSSLNRuHaCYab4kx/umG9muATwBvjTH+CfAVmsnQDcAa8Argln7HK0mSJI4KJfUAAAhpSURBVGl3sliz8VvABPCrIYTPhhB+oNX+ROC7gDGAGOODwDuAz9Ic1XhfjPFc/8OVJEmStBt9H9mIMT53m/ZzNJOQjW0fBT7aj7habu3juTqVt5jyFg8YU6fyFlPe4gFjUm/k9f8wj3HlMSbIZ1zG1Lk8xjU0MfVsgbgkSZKk4eamfpIkSZJ6wmRDkiRJUk+YbEiSJEnqib4vEM+j1m7ltwJPBfYBr4sx/nm2UUEIYQb458AzY4xfzTgcQgjfC7wLuBK4Lcb4wYzjuQr418AB4CrgH8cY/21GcfwEcBw4EWP8cOu1+mWgAJyJMb49BzG9Gng7MAp8Ksb4jqxjarVfDfx34L/GGF+fZTwhhBHg3cD3A48AfzvG2Mg4ph+iWZnvW8BfAq+NMX6zXzFp9/Lat0D++hewj+kgFvuZXcbVas+kr9kupmHqbxzZAGKMF4B3xhhfCXwQeEvGIRFCOAj8EPAHWccCEEIYpbnB4utijD+QdSfQ8iLg4RjjS4A3Aa/PKI4J4JvAbRvabgVuiDF+P/D8EMLhHMT0x8D/AXwv8LrW3jZZxwTwT4Df6HMs0D6efwgsxhhfHmN8ZT8/+C8S09uAH40xvpTmJkrP7nNM2qU89i2Qv/4F7GM6YD9zeXFBdn0NDHl/Y7LREmP8i9ZfDwBns4wlhJAAvwL8FM1sNw9eCjwBOB1C+C8hhNdkHRBwJ3B1COEDND9EfjGLIGKMtRjjB4BVeHSDylXgwRD+//buPUaPqg7j+LcWa6DtKoiCSmyIKxivXL0QgSaoMQ2UxtgnxEuzSmmK0YQmJooQU6tNEaoYpGqLVrmI9DGmXgioESzWEkGNSZUIhYR4ISkBQ1IqtIDWP86pLCu7fXdfZ2e6+3z+2vfMOzO/3bw7T86ZM+/Rt4A5wKlt1lTb/mr7acqI0x5gd9s1STqTMor508msZbR6gHOBd0i6TdL1kmZ3oKa1wPWS1gMPAX+YzJqiP13KFuhsvkAyZkzJmf7qajNrRquJaZQ36WwMI2kRcCKwruVSlgG3297Rch3DDQC/s70AeB9lBKptA8Bs4G7KP8ukXmgP4JXADcBV/O/oSmvqxexG4GM1ENqs5TBgJeWWe1cMACtsnwXsBJa0XA/APGAHcC/l+tTGSGH0oUPZAt3MF0jGTERyprd6upg1MI3yJp2NStL5wGJgse1nWi7nHGCRpC3ACcBNko5vtyTup3wIoczl+2eLtey3ENhuewMlQC9ouR4AbD9GGc1ZSrmlvIAOTFeQdCSwGfii7V+0XQ9wGvBS4EeUuePvlbSm3ZKe8znfS7nF3LZPA5+wfSXwR+CMluuJcehYtkA38wWSMeOSnBmXLmYNTKO8yaJ+gKRTgLuAbZTbyk/Zfk+7VRU1EIa68ACfpK9Q5mE+DayyfVvL9byCMprzAsqXHVxp+wct1HEM8EPKKNNeYCvwHWAN5W/1c9urOlDT4cAbgb/Ut22wfWObNdleUrfNp3zOh9qsB/gy8DXK7f+/AUsnc2RulJp+Sxnx2g08CnzE9qRPTYjx63K2QLfyBZIxB6glOdNHXW1mzWg1MY3yJp2NiIiIiIhoRKZRRUREREREI9LZiIiIiIiIRqSzERERERERjUhnIyIiIiIiGpHORkRERERENOKQtguIGK/6dY3HAU/Upm/avmzY9p8AD9heMWK/y4FbbG8Z0X4kcF1dTGoi9XwcOKWJr9KTdDKwCTgC2Gi7a4sSRURMCcmWiGaksxEHq2W2bx5l24PA35+n/a3Ar56nfU7d1jm2fw8MSlpJqTMiIpqTbIn4P0tnI6YMSUcDvwbmUhan2d++mLLw0auADZKeAHbZPknSTcDbgcMlPVB3+bbt1ZIOBb5KWX10JrDW9jX1mO8ErgBeRhkZ+vEBansDZXXXo2zvkTSDElwfsH2npE3AycAM4B7gg7YfP8Ax9wFzbe+WdDbwSdvzJc0EVlNWCp5FGVn7fN3neGA9cDRwKGWl0DFrj4iYzpItyZboT57ZiCnD9k7bg8CnRrR/v7bfRRm1GrR9Ut12HjAfeKy2D9peXXf9DHCf7ddTQuMSScdKegnl9vPSetzP9VDbPcC9wMLadAbwpO076+vl9dyvAR6nrOA5UUPAbMoqrm8GFko6s267GLjZ9utszwNu6eM8ERFTXrLlv4ZItsQE5M5GHKy+Lmlt/fli25sbOMcC4AhJ59fXs4DXUkai/lwv8gD/6vF4G4EPAwY+VF/v91FJS+o5Xg7c32fdbwPeXV/PAQaBO4BbgS9JGgC+a/u+Ps4TETHVJFvGrjvZEuOWzkYcrC4cY17taPaNs30mMGT7juGNkt4P7B3W1Ov/0feAyyUdA5wLvKke7zRgBXCi7UckXTbGMXoxE7jE9rUjN9jeJOlu4DzgZ5KusL2uz/NFREwVyZbRJVtiQjKNKqaTR4C3ANR5rfv9AxiQ9OoR27YCF0l6UW0/rLZvB06V9GJJxwEX9XJy27so82/XA9tsP1w3HVVre1TSIM/eDu/l95knaRawfFj7VmB5HWFC0gslHVJ/nmv7QdtrgC/w7AhVRERMTLIl2RJjyJ2NmDIkvQv4BuUhvtn1obyVtm+ob1kDXCdpKbCLGg71IbhLgW2SngKuBVYBnwWuBnZIehLYA5xge4ek9cCfKA/iXUWZv9qLjcDtPPeifyuwrB7rIeA3w36nIeBSyoOCMyQtAi6w/UtgJbAZeBjYApxed7saOBbYLmkP8AxwVn3fOkmnA/8GdgIX9lh3RMS0lGxJtkR/ZuzbN9pdvoiIiIiIiInLNKqIiIiIiGhEOhsREREREdGIdDYiIiIiIqIR6WxEREREREQj0tmIiIiIiIhGpLMRERERERGNSGcjIiIiIiIa8R/J4OGRIMrsPwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 估计值\n",
    "homo_x = homo_fitted_result.fittedvalues \n",
    "hetero_x = hetero_fitted_result.fittedvalues \n",
    "\n",
    "# 残差值\n",
    "homo_e = homo_fitted_result.resid\n",
    "hetero_e = hetero_fitted_result.resid\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 6))\n",
    "fig.tight_layout(pad=5)\n",
    "\n",
    "axes[0].scatter(homo_x, homo_e, color='b')\n",
    "axes[0].set_xlabel(\"Fitted values\")\n",
    "axes[0].set_ylabel(\"Residual\")\n",
    "axes[0].set_title(\"Residuals vs. Fits Plot\")\n",
    "\n",
    "axes[1].scatter(hetero_x, hetero_e, color='r')\n",
    "axes[1].set_xlabel(\"Fitted values\")\n",
    "axes[1].set_ylabel(\"Residual\")\n",
    "axes[1].set_title(\"Residuals vs. Fits Plot\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "表1\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>均值(homo)</th>\n",
       "      <th>方差(homo)</th>\n",
       "      <th>均值(hetero)</th>\n",
       "      <th>方差(hetero)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>x=1</th>\n",
       "      <td>0.234120</td>\n",
       "      <td>0.795501</td>\n",
       "      <td>1.481288</td>\n",
       "      <td>28.638050</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=3</th>\n",
       "      <td>0.562194</td>\n",
       "      <td>1.467287</td>\n",
       "      <td>4.780266</td>\n",
       "      <td>93.906343</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=5</th>\n",
       "      <td>0.381242</td>\n",
       "      <td>0.531841</td>\n",
       "      <td>4.371954</td>\n",
       "      <td>53.184129</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=7</th>\n",
       "      <td>0.170099</td>\n",
       "      <td>0.637128</td>\n",
       "      <td>2.948230</td>\n",
       "      <td>91.746470</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x=9</th>\n",
       "      <td>-0.557512</td>\n",
       "      <td>0.810359</td>\n",
       "      <td>-6.479940</td>\n",
       "      <td>158.830340</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     均值(homo)  方差(homo)  均值(hetero)  方差(hetero)\n",
       "x=1  0.234120  0.795501    1.481288   28.638050\n",
       "x=3  0.562194  1.467287    4.780266   93.906343\n",
       "x=5  0.381242  0.531841    4.371954   53.184129\n",
       "x=7  0.170099  0.637128    2.948230   91.746470\n",
       "x=9 -0.557512  0.810359   -6.479940  158.830340"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "homo_nd = homo_e.reshape(5, 10)\n",
    "homo_df = pd.DataFrame(homo_nd) \n",
    "homo_mean = homo_df.apply(np.mean, axis=1)\n",
    "homo_var = homo_df.apply(np.var, axis=1)\n",
    "\n",
    "hetero_nd = hetero_e.reshape(5, 10)\n",
    "hetero_df = pd.DataFrame(hetero_nd) \n",
    "hetero_mean = hetero_df.apply(np.mean, axis=1)\n",
    "hetero_var = hetero_df.apply(np.var, axis=1)\n",
    "\n",
    "print(\"表1\")\n",
    "pd.DataFrame(\n",
    "    data=np.column_stack((homo_mean, homo_var, hetero_mean, hetero_var)),\n",
    "    index=[\"x=1\", \"x=3\", \"x=5\", \"x=7\", \"x=9\"],\n",
    "    columns=[\"均值(homo)\", \"方差(homo)\", \"均值(hetero)\", \"方差(hetero)\"])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
